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Summary 

Delayed division is an iterative method for solving the linear eigenproblem AX = ABX for a limited 
number of small eigenvalues and their corresponding eigenvectors. The distinctive feature of the method is 
the reduction of the problem to an approximate triangular form by systematically dropping quadratic terms 
in the eigenvalue A. The report describes the pivoting strategy in the reduction and the method for preserving 
symmetry in submatrices at each reduction step. Along with the approximate triangular reduction, the report 
extends some techniques used in the method of inverse subspace iteration. Examples are included for problems 
of varying complexity. 

Introduction 

A step in the postbuckling analysis of plates and shells is the computation of several small eigenvalues 
and corresponding eigenvectors of a linear eigenproblem AX = ABX. The matrices A and B that appear in 
the postbuckling problem are symmetric but not necessarily positive definite. One condition assumed in most 
methods for solving the eigenproblem is that either matrix A or matrix B is positive definite, and therefore 
most standard computer subroutines are not applicable for solving the eigenproblems that arise in the solution 
of postbuckling problems. A new algorithm for solving the symmetric eigenproblem has been developed as an 
integral part of postbuckling analysis. 

The new algorithm, called “delayed division,” solves the symmetric problem AX = ABX and is described in 
this report. The contents of the report are confined to the eigenproblem per se, independent of its application to 
postbuckling analysis. This report also describes the general properties that make delayed division a suitable 
choice as a method in the solution of other applied problems as opposed to the choice of the method of 
determinant search or the method of subspace iteration. 

The method of delayed division was first applied to computing roots of lambda matrices (ref. 1). In 
connection with that problem, two computer subroutines were written for solving the problem AX = ABX, 
where A and B are arbitrary square matrices, for only the eigenvalue nearest to zero and its corresponding 
eigenvector. One subroutine searched for only real eigenvalues, whereas the second subroutine computed 
complex eigenvalues. The subroutine for real eigenvalues has been applied by Anderson (ref. 2) to the buckling 
of periodic lattice structures. 

Delayed division is an iterative method. As with other iterative methods for eigenproblems, its rate of 
convergence and its efficiency are problem dependent. This report contains information that pertains to the 
efficient application of the method and describes the following: the basic algorithm for delayed division; the 
logic of a computer subroutine (SPLIT) that was used to implement the algorithm; simple numerical examples 
that illustrate the computer logic; and example problems that were solved by applying the SPLIT subroutine. 

The basic algorithm is an extension of the method of Gaussian elimination for linear algebra problems. 
Once an eigenvalue is known, Gaussian elimination can be applied to the matrix [A — AB] to compute the 
corresponding eigenvector. The reduction of the matrix to triangular form involves a number of divisions. The 
method of delayed division is based on a reduction of the matrix [A — AB] when the eigenvalue is not known. 
The divisions in the reduction are approximated by a truncated power series. The approximate divisions are 
completed after approximate eigenvalues are computed from the truncated triangular form. The completion of 
the reduction in a later step of the computations is the reason for the name delayed division. 

The report describes the logic of the new SPLIT subroutine that was specifically designed for the case 
in which matrices A and B are symmetric and the iteration procedure is directed at computing more than 
one small eigenvalue at a time. As in Gaussian elimination subroutines, the SPLIT subroutine has a forward 
sweep subroutine and a back substitution subroutine that contain most of the arithmetic operations. The 
forward sweep reduces matrix A to exact triangular form and matrix B to approximate triangular form. After 
computing approximate eigenvalues from the triangular forms, the corresponding eigenvectors are computed 
in the back substitution. The approximate eigenvalues identified in the forward sweep can be corrected by use 
of an eigenvalue shift and repetition of the forward sweep and back substitution. 

The basic algorithm of the method of delayed division consist of a forward sweep followed by a back 
substitution or by an eigenvalue shift and can be viewed as an extension of the method of determinant search. 
In addition to the basic algorithm, an optional strategy in the SPLIT subroutine is a modified shift that 
corrects more than one eigenvalue and corresponding eigenvector at a time. This option is similar in scope to 
the method of inverse subspace iteration. 



Discarding the assumption that either matrix A or matrix B is known to be positive definite introduces the 
need for a pivoting strategy in the solution options. The computer logic for the pivoting strategy and for the 
modified shift is described in detail in the section entitled “Delayed Division.” Each major step in the algorithm 
for delayed division is contained in a SPLIT subroutine. The subroutines are discussed in the order that they 
are called by the main subroutine starting with the forward sweep procedure. The numerical operations in 
each step are illustrated with simple numerical examples. 

After the description of the computer logic in the SPLIT subroutine, results are presented from four examples 
that apply delayed division to the symmetric eigenproblem AX = ABX. The first example, the finite difference 
equations for a column buckling problem, provides a measure of the rate of convergence of the iteration as 
a function of both the number of intervals in the difference equations and the computed number of critical 
loads. The second example is a textbook problem that demonstrates the numerical stability in the eigenvector 
computations; the stability results from the pivoting strategy in the forward sweep subroutine. The third 
example is the shear buckling of simply supported plates, an example which shows the user can often modify 
problems to speed up convergence. The fourth example is an application of delayed division to lambda matrices 
derived by the method of Lagrangian multipliers, a practical application that combines the power and the 
versatility of the method of Langrangian multipliers with the rapid convergence of delayed division. After the 
examples, the report concludes with a discussion of the similarities and differences between delayed division 
and other methods for solving the symmetric eigenproblem. 


Delayed Division 


Basic Algorithm for General Case 

The method of delayed division is a simple extension of the fundamental idea in Gaussian elimination. This 
section reviews the algorithm in reference 1 as applied to the problem AX = ABX, where A and B are arbitrary 
square matrices. This review helps to define the requirements for extending the computer implementation of 
the method to the symmetric eigenproblem. 

The basic algorithm in delayed division follows from the steps in the computation of the eigenvector X by 
the method of Gaussian elimination applied to the equation (A— AB)X = 0. The elimination procedure requires 
division by terms of the form (a + A b), where a and b are constants. Division by (a + A b) is approximated for 
small eigenvalues by multiplication of the two leading terms in the power series for the reciprocal of (a + A6) 


_1 _1_A6 (A b/a) 2 

(a + A6) a a 2 + (a + A b) 


( 1 ) 


1 1 _ Xb 

(a + Xb) a a 2 



(2) 


The two coefficients (1/a) and (—b/a 2 ) are stored; thus, division by (a + Xb) when A is unknown can be delayed 
until a value for A is assigned. 

The approximation in equation (2) is used in reference 1 in a forward sweep procedure to produce an 
approximate upper triangular matrix. The matrices A and B are read into auxiliary matrices C and D such 
that C = A and D = — B. The forward sweep operates on these auxiliary matrices. A typical operation 
consists of solving for X/. in an equation of the form 


{ c ik T Xd ik )X k — ^ ^ (c{j + Xdij )Xj 

The exact solution of equation (3a) for ~K k in terms of the Xy is 


(3a) 




-E 

J^k 


( c ij 

( c ik “t~ Xdi k ) 


(3b) 
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The approximate solution of equation (3a) by use of the right-hand side of equation (2) as the divisor and 
without quadratic terms in A in the multiplications is 

X fc--E(4+H) x i ( 3c ) 


where 



j _ ^*3 
13 c ik 


c ijdik 

~Aj2 

c ik 


The terms c( • and dG are stored over c^j and d t j . Equation (3c) is used to eliminate X& from the remaining 
equations; the usual algebraic rules of Gaussian elimination are applied except that quadratic terms in A 
are dropped as they appear. An elementary numerical example will show the simplicity of delayed division 
and, at the same time, will reveal details of numerical analysis that must be kept in mind in the computer 
implementation. The numerical example is a textbook problem from reference 3, wherein it is used to illustrate 
convergence of the power iteration for the standard problem AX = AX. The problem is to find nontrivial 
solutions of the following homogeneous equations: 


(3 - A)Xi + 2X 2 + 1X 3 = 0 (4a) 

2Xi + (2 - A)X 2 + IX3 = 0 (4b) 

lXi + 1X 2 + (1 — A)X 3 = 0 (4c) 

The recursive use of equation (2) follows along the lines of Gaussian elimination, with quadratic terms in 
A dropped at each step in the elimination process. Division of equation (4a) by the coefficient (3 — A) is 
approximated from multiplication by (1/3) + (A/9) to give 

lXi + [(2/3) + (2A/9)]X 2 + [(1/3) + (A/9)]X 3 = 0 

The next step is elimination of Xi in equations (4b) and (4c) by multiplying the above equation by the 
appropriate coefficient of Xi, dropping quadratic terms in A, and subtracting. The result is 


lXi + [(2/3) + (2A/9)]X 2 + [(1/3) + (A/9)]X 3 = 0 (5a) 

0X x + [(2/3) - (13A/9)]X 2 + [(1/3) - (2A/9)]X 3 = 0 (5b) 

OXi + [(1/3) - (2A/9)]X 2 + [(2/3) - (10A/9)]X 3 = 0 (5c) 

Multiplication by (3/2) + (13 A/4) is the first step in eliminating X 2 from equation (5c). The resulting 
equations are 

lXi + [(2/3) + (2A/9)]X 2 + [(1/3) + (A/9)]X 3 = 0 (6a) 

OXi + 1X 2 + [(1/2) + (3A/4)]X 3 = 0 (6b) 

OXi + 0X 2 + [(1/2) + (5A/4)]X 3 = 0 (6c) 


For a nontrivial solution, the coefficient of X 3 in equation (6c) must become zero. The coefficient is zero when 


A = Ai = 2/5 = 0.40 (7) 

However, the triangulated equations (6) are not exactly equivalent to equations (4), since quadratic terms in 
A have been dropped at each step. The algorithm, through use of equation (2), continues the iteration by 
application of the following eigenvalue shift fi in equations (4): 

A = 0.40 + n (8) 


3 



After the shift, equations (4) are 


(2.6 - /z)X i + 2X 2 + 1X 3 = 0 (9a) 

2Xi + (1.6 - //)X 2 + IX3-O (9b) 

lXi + 1X 2 + (0.6 - //)X 3 = 0 (9c) 

The forward sweep through equations (9) shows the need for pivoting before delayed division can be expected 
to converge in the general case. After the first elimination step for equations (9), the resulting equations are 

lXi + (0.7692 + 0.2959//)X 2 + (0.3846 + 0.1479//)X 3 - 0 (10a) 

OXi + (0.0615 - 1.592//)X 2 + (0.2308 - 0.2959//)X 3 - 0 (10b) 

OXi + (0.2308 + 0.2959//)X 2 + (0.2154 - 1.148//)X 3 = 0 (10c) 


For a given value of //, the approximation of equation (2) is more accurate for the reciprocal of (0.2308— 0.2959/t) 
than for the reciprocal of (0.0615 — 1.592//), since the absolute value of b/a is smaller for the first expression. 

Therefore, a proper choice is to solve for X 3 in equation (10b) and eliminate X 3 in equation (10c). After the 

complete forward sweep, the equtions are 

lXi + (0.7692 + 0.2959 //)X 2 + (0.3846 + 0.1479//)X 3 - 0 (11a) 

OXi + (0.2667 - 6.556 //)X 2 + 1X 3 = 0 (lib) 

OXi + (0.1733 + 1.422//)X 2 + 0X 3 = 0 (11c) 

For a nontrivial solution of equations (11), 

H = -0.121875 (12) 

A = 0.40 + // = 0.278125 (13) 

For this value of /t, the approximate reciprocal of (0.2308 — 0.2959//) calculated from equation (2) is correct 
to within 3 percent, whereas the approximation for (0.0615 — 1.592//) has the wrong algebraic sign, since 
| b/z / a | > 1. Therefore, the use of the approximate reciprocal from equation (2) as part of a numerically stable 
forward sweep algorithm requires permuting the order of elimination of the variables, or “pivoting for size.” 
(See ref. 4.) The algorithms used in the subroutines described in reference 1 use “maximal pivots” determined 
by selecting the maximum absolute value of Cy from the coefficients (Cy + Ady ) that remain in the equations 
of lower order after each elimination step. Maximal pivoting on the off-diagonal terms destroys the symmetry 
of the reduced system of equations of lower order when A and B are symmetric matrices. Preservation of 
symmetry to decrease the total number of calculations is the motivation for using a different pivoting strategy 
for symmetric problems. 

Before another pivoting strategy is presented, the convergence of the iteration calculated with the 
approximate reciprocal and an eigenvalue shift is summarized below. 

aJ 1} - 0.40 A ( j 4) = 0.3079788 

A? 5 = 0.278125 aS 5) = 0.3079785 

A ( j 3) = 0.3085412 Exact Ai = 0.30797853 

Note that the eigenvector associated with A can be computed by back substitution in the triangulated equations 
when //, the correction on A, is small. 

Approximate 2x2 Matrix Inverse 

The pivoting strategy in the algorithm in reference 1 destroys the symmetry of the reduced system of 
equations of lower order when an off-diagonal term is the pivot element in an elimination step. This section 
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introduces another pivoting strategy that preserves the symmetry of the original problem. The new pivoting 
strategy is based on eliminating two variables at a time while retaining a systematic procedure for dropping 
quadratic terms in A as they appear. Eliminating two variables at a time involves computation and storage of 
an approximate inverse of a symmetric 2x2 matrix. If the matrix is written as 


the exact inverse is 


(C + AD) 


(('li T A (i 7 j ) ( c^j -f - Xdjj j 

(cjj + Xd^j ) ( Cjj + Xdjj) 


(C + AD)" 1 = 


( c jj + Xdjj) /A 
“h Xd^j)/ A 


(c ij Xd^j ) f A 
( c n "f Xd#)/A 


which is also symmetric. 


In equation (15), the divisor A is the determinant of (C + AD): 


(14) 


(15) 


A — cl bX 4” cA 2 


(16) 


where a = ( c^cyy - c?-)> 
in delayed division is 


b = 


{ c iidjj ~t~ Cjjda 

(C + AD)- 1 


2Cydjj), and c ( d^djj d^y). The 

( c 'ii + Xd' i{ ) (c'y + Ad'y) 

( C ij + Xd ij) [c’jj + Xd'jj) 


approximate inverse used 


(17) 


where 


J C 33 


j _ c ij 
a 



d'ii 

1 

tTI 

1 

b c jj 

a 

a 2 

<•; 

4 

1 

1 

be, 

a 

i : 

d‘ 

d 'n 

d ii 


a 

a 2 


The approximation in equation (17) follows from dropping the quadratic term in A in equation (16) and then 
applying the approximation of equation (2) to the division by A in equation (15). The approximation preserves 
the symmetry of the inverse. The approximate inverse in equation (17) is the unique feature of delayed division 
for the symmetric problem AX = ABX. Before the approximation in equation (17) is applied to the general 
elimination procedure, it will be used on the simple numerical example in equations (4). The leading 2x2 
coefficient matrix in equations (4) (e.g., A 2 ) corresponds to equation (14) as follows: 


" (3 - A) 2 ‘ 

-1 

(2 — A)/A 

-2/A ‘ 

2 (2 — A) 


-2/A 

(3 — A)/ A . 


where A = (3 — A) (2 — A) — 4 = 2 — 5A + A 2 . The approximate inverse of A 2 is 


(18) 



‘ (1 + 2A) — 1 — (5A/2) ' 

. — 1 — (5A/2) (3/2) + (13A/4) . 


(19) 


Equation (19) follows from equation (18) with (1/A) approximated by (1/2) + (5A/4) and the quadratic 
terms in A dropped. After multiplication of equations (4a) and (4b) by Ag \ elimination of X] and X 2 from 
equation (4c), and elimination of quadratic terms in A at each step, the resulting equations are 

lXi + 0X 2 + (0 - A/2)X 3 = 0 (20a) 

OXj + 1X 2 + [( 1/2) + (3A/4)]X 3 = 0 (20b) 

0X l + 0X 2 + [(1/2) - (5A/4)]X 3 = 0 (20c) 

Equations (20) and equations (6) are equivalent; they have the same nontrivial solution. The order in which 
quadratic terms are dropped does not affect the final result of the forward sweep. 
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At the next iteration step for the numerical example, the pivoting problem that arises in solving equations (9) 
is resolved by first solving the equations 


(2.6 - n)X i + 1X 3 - -2X 2 (21a) 

1X X + (0.6 - /z)X 3 = -1X 2 (21b) 

for Xi and X 3 in terms of X 2 . The choice of order of elimination is based on the determinant for the coefficients 
on the left side of equations (21). That determinant is A = 0.56 — 3.2A/i + A/i 2 . whereas the determinant for 

the leading 2x2 coefficient matrix of equations (9) is A = (2.6 — A^u)(1.6 — A/jl) — 4 = 0.16 — 4.2A/Z + A/z 2 . 

The constant term 0.56 is larger than 0.16, so that Xi and X 3 are eliminated first with equations (21). Again, 
the final triangular equations at the end of the forward sweep are equivalent. The symmetry of the reduced set 
of equations at each step of the forward sweep is not readily apparent for the simple 3x3 matrix here. The 
preservation of symmetry is illustrated in the next subsection. 


Forward Sweep 

In this section, three numerical operations from the forward sweep procedure for the symmetric problem 
are outlined in more detail. The first operation is the use of 2 x 2 inverse matrices in the elimination procedure. 
Symmetry in the reduced equations is preserved for the general case. The second operation is the pivoting 
strategy, that is, the choice of the order in which variables are eliminated in the SPLIT subroutine. By knowing 
the pivoting strategy, the subroutine user can often scale applied problems to permute the order that variables 
are eliminated. The effect of scaling problems with equivalent matrices is discussed briefly. The third operation 
is inclusion of quadratic terms in the determinant of the final 2x2 matrix at the end of the sweep. The third 
operation gives two approximate small eigenvalues rather than one, with no additional approximations. 

The first operation considered herein is the elimination of two variables at a time in the forward sweep 
procedure. As in Gaussian elimination for linear algebraic equations, the A matrix is replaced by a sequence 
of equivalent matrices. In the computer algorithm, only the current equivalent matrix is stored. The notation 
adopted here is to let C = A and D = — B at the start of the forward sweep. Before each elimination step, 
the equations are of the form 

(C + AD)X = 0 (22) 

After the elimination step, a prime is used to distinquish the new set of equations: 

(C' + AD')X = 0 (23) 


The terminology “reduced set of equations” is usually applied to the new set of equations (eqs. (23)). 
However, in this report “reduced equations” are the equations in variables that remain to be “eliminated” 
in equations (23). To condense the notation somewhat, let 


9ij — c ij T \d{j 


/ (* = 1 , 2 , •••>») 

\ (j = 1, 2, ..., n) 


where and d{j are elements of the n x n matrices C and D. 

At the first elimination step, equations (22) are symmetric so that gji = g^j . This symmetry is preserved 
in the n — 2 reduced equations after the first step, g^- = g'-. The elimination process is then repeated for the 
reduced set. The symmetry of the reduced equations is apparent when the elimination step is written in detail. 
The equations, before the first elimination step, are of the form 


n 

= 0 (* = l,2,...,n) 

3 = 1 


The operation considered is the elimination of X/ and X j where I and J are selected so that 


CiiCjj-C 2 


IJ 


— c a c jj c ij 


f (*= 1, 2, ..., n) 
\(j = 1,2,..., n) 


( 24 ) 
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In other words, I and J are selected to give the maximum absolute value for the symmetric 2x2 submatrices 
that appear in the coefficient matrix for the equations. The /th and Jth equations of equations (24) are written 


9lI x I + 9IJ X J 


E 9lj x j 

m,J 


(25a) 


9ij x I + 9 jj x j = ~ E 9jj x j 


(25b) 


The exact inverse matrix of coefficients on the left side of equations (25) is defined in a manner similar to 
equation (15) as 

V/ 1 9jJ 


with 


where, as in equation (16), 


GJj = 


9lJ 9JJ J 


9 II = ( C JJ + Mjj)/A 
9jj = -( C IJ + Mu )/ A 
9jj ~ ( c n + Mu)/ A 

A = a + bX 4“ cA 2 


(26) 

(27a) 

(27b) 

(27c) 


and a = {c n djj -cjj), b = {cndjj + cjjd u -2c/jd/j), and c = {d n djj -djj). Multiplying equations (25) 
by the inverse matrix of equations (26) yields 


xi 


= ~ E & 3 X 3 


3*I,J 


X J = ~ E 9jj> 


(28a) 

(28b) 


where 


}&,J 

9ij = 9 u 9lj + 9ij9jj 
9jj = 9u9lj + 9jj9jj 

When i ^ I or i ^ J in equations (24), the typical equation with the row subscript i is 

gu x i + guxj + E gijXj = o (*V ^ J) 


(29) 


The variables xj and x j are eliminated from equations (29) by substitution of the right side of equations (28) 
to give 

E 9ij x j — 0 (* I, * 7" J) (30) 


where 


or 


9ij ~ 9ij 9il9lj 9ij9jj 


g'ij = 9xj ~ 9il(9ll9lj + 9lj9Jj ) - 9ij(9lj9lj + 9jj9Jj ) 
In equation (32), interchanging subscripts shows directly that 


(31) 

(32) 


J — J 

9ji 9ij 


( 33 ) 
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Therefore, the reduced equations (30) are symmetric. In the computer subroutine, the reduced-equation 
coefficients are truncated to linear terms in A and then the c'- and d'- terms are stored over the Cjj and djj 
coefficients. Each term in equation (32) is a product or quotient of elementary linear factors of the form (c+Ad) 
so that the truncation to linear terms in A can follow any order. The arbitrary order of truncation makes it 
convenient to truncate equations (27) first. The coefficients in the term 

<*> 


are stored in the original locations for c/j and djj and the coefficients of gj j and gjj are stored in a similar 
manner. 

The purpose of preserving symmetry in the reduced equations is to decrease the number of operations in 
the numerical solution. Only half the off-diagonal terms are computed at each elimination step. 

The second operation in the forward sweep that is considered in more detail herein is the pivoting 
strategy. The variables xj and xj to be eliminated at each step are arbitrarily chosen so that the coefficient 
a = (c/jc/j-cfj) has a maximum absolute value of all possible (c&cyy — c|-). The accuracy of the approximate 

truncation to linear terms of the inverse matrix Gjj defined by equations (26) and (27) depends on an 
application of equation (7). The reciprocal for A, defined in equation (16), is 


A 


1 

a 


1 — A 



.2 [(&/«) + (*c/a)] 2 
+ A - 


(35) 


For the leading term in equation (35) to 
satisfy the inequality 


The accuracy of the truncated expression 


approximate the reciprocal of A, the term (A b/a) + (A 2 c/a) must 



depends on A, b, and c as well as the constant a. However, the algorithm programmed in the SPLIT subroutine 
only tests for the magnitude of a. Since A is not known, there is no control over A during the forward sweep. 
Before the sweep procedure starts, an eigenvalue shift can be used to speed up convergence. 

The ratios (b/a) and (c/a) appearing in inequality (36) could be used as a basis for a test to select pivot 
variables xj and x j. Such a test would depend on matrix elements djj in addition to Cjj and would require 
more computer operations per elimination step, so that this additional testing is neither explored further herein 
nor adopted in the SPLIT subroutine. However, the subroutine user can exert some control over the sweep 
process by scaling the original problem. One means of scaling the symmetric problem (A — AB)X = 0 is to 
change X to RY and multiply the matrix equation by H T , the transpose of R. The original symmetric problem 
becomes a symmetric equivalent problem 


(R r AR)Y = A(R r BR)Y (38) 

with the same eigenvalues. The matrix C at the start of the forward sweep is R T AR. The order of elimination 
of variables for the matrix C in the subroutine can obviously be different from that of the original matrix A. 

Fortunately, most mechanics problems tend to scale naturally. Scaling is illustrated in this report for the 
problem which applies the method of Lagrangian multipliers to plate buckling. 

The pivoting strategy adopted herein is numerically stable and the elimination procedure is exact when A is 
zero. The operations in which matrix A is factored contain no approximations. If the interchange of row and 
column matrix elements is included in the forward sweep procedure, the final matrix C is triangular and the 
product of the (c//Cj j — c|j) terms generated by the sweep are equal to the determinant of A. Interchanges 
of elements are not included in the forward sweep procedure in the SPLIT subroutine, but the value of the 
determinant is independent of interchanges and is computed as a product after the sweep. 
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The third operation in the forward sweep occurs at the last elimination of the sweep, when the reduced 
set of equations contains two equations. In this case, the eigenvalues of the determinant of the last 2x2 
reduced matrix are computed exactly. The computation gives two eigenvalues that are approximations of 
small eigenvalues of the original problem. 

This last operation is not illustrated in the elementary numerical example considered previously; the 3x3 
matrix of that example is too small to show the practical significance. The following augmented problem is a 
better example: 

(3 - A)Xi + 2X 2 + 1 X 3 + OX 4 = 0 (39a) 

2Xi + (2 - A)X 2 + IX 3 + OX 4 = 0 (39b) 

1 X 4 + 1 X 2 + (1 - A)X 3 + OX 4 = 0 (39c) 

OXi + 0X 2 + OX 3 + (4 - A)X 4 = 0 (39d) 


The numerical example defined by equations (39) has eigenvalue A 4 = 4 plus the eigenvalues of the first 


numerical example. 

The forward sweep procedure starts with 7 = 1 and J = 4. To eliminate Xi 
(eqs. (26)) is 



r_[4 z A) 

12-7A+A 2 

(0 + 0A) 


(0 + 0A) 
(3— A) 

12-7A+A 2 J 


■(1/3) + (A/9) 
. (0 + 0A) 


(0 + 
( 1 / 4 ) + 


and X 4 , the inverse matrix 


0A) 

(A/16). 


(40) 


After the elimination of Xi and X 4 in the reduced equations, equations (39) become 


lXi + [(2/3) + (2A/9)]X 2 + [(1/3) + (A/9)]X 3 + 0X 4 = 0 (41a) 

OXi + [(2/3) - (13A/9)]X 2 + [(1/3) - (2A/9)]X 3 + 0X 4 - 0 (41b) 

OXi + [(1/3) - (2A/9)]X 2 + [(2/3) - (10A/9)]X 3 + 0X 4 = 0 (41c) 

OXi + 0X 2 + OX 3 + IX 4 = 0 (41d) 

The reduced equations in X 2 and X 3 are equations (41b) and (41c), which are symmetric. For a nontrivial 
solution, the determinant A of the coefficients must vanish: 

A = [(2/3) - (13A/9)][(2/3) - (10A/9)] - [(1/3) - (2A/9 )] 2 = 0 (42) 


The roots of equation (42) are 

Xi = (7 — \/7) /14 = 0.3110 (43a) 

A 2 = (7 + Vl) / 14 = 0.6890 (43b) 

The exact eigenvalues of the given problem (eqs. (39)) are Aj = 0.30797853, A 2 = 0.64310413, A 3 = 5.04891734, 

and A 4 = 4. The roots of equation (42) approximate the first two eigenvalues of the problem. In particular, 

the approximation for Ai = 0.3110 is much better than the value Ai = 0.40 obtained from equations ( 6 ) for the 
first numerical example. Equations ( 6 ) follow from dropping linear terms in A in the forward sweep until only 
one reduced equation remains. For the second numerical example, there are two small eigenvalues and two 
larger ones and delayed division provides a good approximation to the first two since the approximate inverse 
matrix in equations (40) is accurate for A = Ai and A = A 2 . 

In the general case, the determinant (eq. (16)) of the final 2 x 2 set of equations is of the form 


A — a + 6 A + cA 2 


with roots 


A = 


—6 ± ( 6 2 — 4 ac ) 1 / 2 
2 c 


(44) 
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As noted in the discussion on the pivoting strategy for selecting I and J in equations (25), the constant a is 
an exact factor of the determinant of the matrix A. Therefore, if a = 0, the root A = 0 of equation (16) is an 
exact eigenvalue of (A — AB)X = 0. When a^O and the roots of equation (44) are real and small, the roots 
provide approximations to two real eigenvalues for (A — AB)X = 0. 

When the discriminant ( 6 2 — 4<xc) is negative, the roots of equation (44) are complex numbers and the 
SPLIT subroutine returns an error message. This raises the question of whether the original problem has 
complex eigenvalues. An eigenvalue shift to the real part of the complex conjugate eigenvalues allows a more 
precise check for complex eigenvalues, or the complex algorithm of reference 1 can be applied. 

Computation of the roots of the determinant of the last pair of reduced equations with equation (44) 
completes the forward sweep procedure. Use of approximate 2x2 inverse matrices reduces the number of 
computations needed in the sweep because of symmetry. Complete pivoting provides a stable computation of 
the determinant of the matrix A . The pivoting also provides a systematic method of dropping quadratic terms 
in A, although this strategy may not be optimum. 

The two roots of equation (44), Ai and A 2 , are small approximate roots of the full eigenproblem. 
Approximate eigenvectors can be computed by back substitution similar to back substitution in Gaussian 
elimination. 


Back Substitution 

This section contains the back substitution procedure for the second numerical example and summarizes 
the matrix notation for the forward sweep and for the back substitution. Two approximate eigenvectors are 
computed by back substitution for the second numerical example. The general case allows for m approximate 
eigenvectors to be computed. The m eigenvectors correspond to m approximate eigenvalues identified in the 
forward sweep. The method of identification is outlined in this section. 

Some additional notation is introduced in this section. The approximate eigenvector corresponding to the 
approximate eigenvalue is denoted as V k . The components of are v kj : 

Vfc = (t>fci,t>*2,~,t;jy) (i = l,2,...,n) (45) 

The notation is reserved for the exact eigenvectors that satisfy the original problem: 

AU fc = AfcBUfc (46) 


The second numerical example is defined by equations (39). At the end of the forward sweep procedure, the 
approximate equations are equations (41). Equations (41b) and (41c) are decoupled equations in X 2 and X 3 . 
For a nontrivial solution, the two equations must be linearly dependent. The linearly dependent equations for 
the components of the eigenvector corresponding to Ai = 0.3110 are written as 


-[(2/3) - (13A!/9)j "[(1/3) — (2A]_/9)] 

13 [(1/3) - (2A!/9)] Ul2 [(2/3) - (10Ai/9)] 12 

V13 = — 0. 8229^12 


(47) 


The arbitrary constant in the eigenvector is fixed by setting V 12 = 1- The remaining two equations (41) 
determine Un and v\$ as follows: 

ln n + 0 ui 4 = - {[(2/3) + (2A 1 /9)]v 12 + [(1/3) + (Ai/9)]« 13 } (48a) 


0vn + Inn = — (0vi2 + OU 13 ) (48b) 

Substitution of Ai,vi 2 , and U 13 into equations (48) leads directly to v\i — —0.43305 and = 0. 

The first eigenvector of equations (41) corresponding to Ai = 0.3110 is, therefore, 
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vf = (-0.4331, 1.0, -0.8229, 0) 


(49) 



The second eigenvector for A 2 = 0.6890 is 


Vj = (0.8596, 0.5486, 1.0, 0) (50) 

The computations for the back substitution in the method of delayed division are straightforward. The implicit 
division that is approximated in the forward sweep is completed by explicit substitution of A& in the back 
substitution. 

The algorithm for improving the approximation of m — 2 eigenvalues and eigenvectors for the second 
numerical example by the use of Vi and V 2 is completed in the next section. Improving the approximation 
for the general case of m eigenvalues, where m > 2, is discussed further here. 

The procedure for selecting m > 2 approximate eigenvalues is part of the forward sweep procedure. For 
each set of equations in the form of equations (25), two eigenvalues are computed from the matrix of coefficients 
on the left-hand side. The matrix of coefficients is 


G u = 


91 1 

9IJ 


.91 J 

9JJ . 



( cji + Ad//) 
(c/j + A djj) 


(c/j + A d/j) 

( C JJ + Mjj) 


(51) 


The determinant A of Gj j has two eigenvalues determined by the roots of 

A = (c// + A d//)(cjj + A djj) - ( cjj + A d/j) 2 = 0 (52) 


The eigenvalues in equations of the form of equation (52) are numbered in reverse order of their computation. 
That is, on the rth elimination step the two roots of the 2x2 matrix determinant are stored as A n +2_2r and 
A n -f-i_ 2 r- At the end of the forward sweep, the m approximate eigenvalues are selected from the A, t with the 
m lowest subscripts. 

The first step in improving the m eigenvalues is the computation of approximate eigenvectors from back 
substitution. The order of the operations in the back substitution is the reverse of the elimination order in the 
forward sweep. In the computer subroutine, a permutation vector P with components is used to record 
the order of operations. When X/ and Xj are the rth pair of unknowns eliminated in the forward sweep, 
two components of the permutation vector are defined as i^r-l — I and P^ r — J. The permutation vector 
determined by the pivoting strategy for the second numerical example (eqs. (39)) is 


P 7 " = (1,4, 2,3) (53) 

The permutation vector is used in the back substitution, but substitution is simpler for the case in which 
the variables are eliminated in their natural order of P \ = i. For this special case, let k = n + 1 — 2r. Then, for 
a nontrivial solution of equations (25) when A = A^, the components of the approximate eigenvector V/. obey 
the relations 


v kl 


n 

~{{ C IJ + A k d Ij) v kJ + E 9lj v kJ ] 
/=2r+l 

(c// + A kdjj) 

n 

+ X k d Jj) v kJ + E 9jj v kj\ 

j=2r+l 


(C/J + A k du) 


(54) 


where I = 2r — 1 and J = 2r. The information necessary to compute v^j for j > 2 r is contained in quadratic 
and higher terms in A that have been dropped. Therefore, the components v^j for j > 2 r are set equal to 
zero at this stage of the iteration. Either v^j or v^j is set equal to one and equations (54) define v^j or v^j 
with absolute value less than one. The remaining v^j for j < 2r — 2 follow by direct back substitution in the 
triangular equations formed during the forward sweep. 

At the end of the forward sweep and the back substitution, m approximate eigenvalues and eigenvectors have 
been identified. The accuracy of the approximation is problem dependent; however, the smaller eigenvalues in 
absolute value are more accurate. One optional solution strategy to improve the approximations is to shift to 
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a given approximate eigenvalue and repeat the forward sweep and the back substitution. A second optional 
solution is to assume that the m eigenvectors are useful approximations that can be improved without repeating 
all the operations in the forward sweep procedure. 

The two solutions are indicated in a schematic diagram in figure 1. The iteration in the first solution makes 
an outer loop that repeats the forward sweep procedure. The optional second solution contains an additional 
inner loop that applies a modified eigenvalue shift. The steps in the inner loop are discussed next. 

Equivalence Transformation 

The forward sweep and back substitution procedures identify m approximate small eigenvalues and their 
associated eigenvectors. This section outlines an equivalence transformation that improves the approximations. 
In addition to improving the approximate eigenvalues, the equivalence transformation allows a modified forward 
sweep on the transformed problem. The modified forward sweep and the back substitution require fewer 
operations than the full forward sweep and the back substitution. The modified forward sweep and back 
substitution procedures complete the first iteration cycle. The modified sweep procedures are described in the 
next section. 

The equivalence transformation for the original problem 


AX = ABX 

(55) 

consists of a change of variables 


X = QY 

(56) 

Followed by multiplication by Q, T . The transformed problem is thus 


(Q t AQ)Y - A(Q t BQ)Y 

(57) 


The original problem and the transformed problem are called equivalent because they have the same eigenvalues. 
A completely general equivalence tansformation allows premultiplication by any conformable matrix. The 
choice of Q for premultiplication makes the transposed problem symmetric, since A = A T and B = B T . 

The reason for transforming the problem is to make it simpler in some sense. The canonical form of the 
change of variables is the special case of Q = U, where the columns of U are the normalized eigenvectors U*. 
of the original problem. The normalized eigenvectors obey the orthogonality relation 

U^Bl Jj = 6 kj (58) 

where 6^ = 0 if k ^ j and 8^ = 1. The choice of Q = U makes the transformed problem (eqs. (57)) 

AY = AIY (59) 

The matrix A = U 7 'AU is a diagonal matrix with elements equal to the eigenvalues A^., since the normalized 
eigenvectors are scaled so that U 7 BU — I. 

The canonical (eqs. (59)) is the end result when all n eigenvalues and eigenvectors of equations (55) have been 
determined. In the method of delayed division, m < n eigenvectors are computed at one time. Transformed 
equations with m exact solutions are examined next to provide insight into the procedure for improving m 
approximate solutions. 

Consider the transformation when the last m columns of Q are eigenvectors tfy and the remaining elements 
are zeroes except for the remaining diagonal elements, which are unity. In partitioned form, Q is written 


Q = 


I n—m,n—m j 


0 m,n—m 


u 


n,m 


(60) 
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Input matrices A and B. 



Figure 1. Schematic diagram of procedure for solving symmetric eigenvalue problem AX = ABX by delayed 
division. 
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where the subscripts on the submatrices denote the number of rows and columns. It is useful here to consider 
the computational steps of the transformation in detail since the computer algorithm for the equivalence 
transformation follows the same steps. The matrix products (AQ) and (BQ) have the partitioned forms 


and 


(AQ) = 


A 


71,71 — m 


! (AU) 


7i, m 


(61a) 


(BQ) = 


B 


71,71—771 


i (BU) n , m 
I 


(61b) 


The notation A n ,n-m and B n , n -m is used to indicate that the first (n — m) columns of the matrices A and 
B are unchanged by the multiplication. Each column k of (AU) n>m is of the form 


AU fc = A fc BU fc 


(62) 


If U fc is in column (n + 1 — k) of Q, (AU^) is in the same column of (AQ). Multiplication of equations (61) 
by Q 7 completes the transformation: 


Q r AQ = 


Q t BQ = 


A n—m,n—m 

(AU) n - m,m 

.(AU )l trl -m 

A m,m 

B n—m,n—m 

(BU) n - m ,m 


(BU) 


T 

m,n—m 


I 

I 


I 


m,m 


(63a) 


(63b) 


The first subscript (n - m) in the notation A n -m,n-m and (AU) n _ m ,m is a reminder that premultiplication 
of a matrix by Q 7 does not alter the first (n — m) rows of the matrix. The m x m identity matrix I m ,m is 
based on the assumption that the eigenvectors have been normalized so that 


U^BU fc = 1 


(64) 


For the normalized eigenvectors, the elements of the diagonal matrix A m ,m are the Rayleigh quotients 


Ufc’AUfc 

k UjTBUfc 


(65) 


The presence of the Rayleigh quotients in the transformed equations suggests that the m approximate 
eigenvalues can be improved by an equivalence tranformation. The orthogonality relation (eq. (58)) can be 
used to improve the m approximate eigenvectors with a Gram-Schmidt procedure. The subscripts on the 
m approximate solutions are determined in a definite order by the pivoting strategy in the forward sweep 
procedure; the ordered subscripts provide a definite order for the Gram-Schmidt procedure. 

The ordered subscripts allow the matrix Q to be formed column by column, and at the same time the 
matrices (QA), (QB), (Q 7 'AQ), and (Q r BQ) are formed. The formation of Q can be summarized as follows: 

1. The m column vectors V' are stored in columns j = P n+ i_fc (for k = 1,2, ...,m), where P is the 
permutation vector determined by the pivoting strategy in the forward sweep. 

2. The (n — m) column vectors of Q have zeroes as elements, except they have ones on the diagonal. 

3. The m vectors Y' are defined recursively by a Gram-Schmidt orthogonalization procedure. The procedure 
is 

= Vi (66a) 
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V^ = V 2 + ?12 V i 

(66b) 

m— 1 


v; = v m + E 9kmK 

(66c) 

k = 1 


The coefficients g^j are computed from 


-(V t ) T (BVy) 
9% (vj.) T (Bvy 

(67) 

The denominator in equations (67) is not unity in general, so that 


(V;) T (BV' ) = d kk 6 kj 

(68) 


The V 7 are orthogonal to the (BVy) but are not normalized. 

The approximate eigenvectors A are updated with the Rayleigh quotients for V(, as follows: 


(n) r (Bv' fc ) 


(69) 


In the computer code the primed quantitites A(, and V(. are stored over the corresponding unprimed variables. 
The numerator and denominator of equations (69) appear as diagonal elements of Q r AQ and Q r BQ. The 
subscript of the diagonal element is determined from the element P n+ i_fc of the permutation vector. 

The equivalence transformation for the second numerical example (eqs. (39)) starts with the eigenvectors 
Vi and V 2 from equations (49) and (50). The permutation vector P for the problem is listed in equation (53). 
For this problem, n = 4 and m — 2; therefore, V( goes in column P(n) = 3 and Q and V 2 goes in column 
P(3) = 2. The change in variables is 


■xr 


'1 -0.8369 -0.4330 0‘ 


vr 

X 2 


0 0.4960 1.000 0 


y 2 

X 3 


0 1.0432 -0.8229 0 


y 3 

.X 4 . 


.00 0 1. 


.y 4 . 


The constant g \ 2 from equations (67) for the second example is 


-0.09796 

1.8646 


-0.0525 


The transformed problem [Q T (A — AB)Q]Y = 0 is 


(3 - A) 

(-0.4752 + 0.8369 A) 

(-0.1220 + 0.4331A) 

0 - 


-Yi- 


- 0 - 

(-0.4752 + 0.8369A) 

(1.30987 - 2.035A) 

(-0.0106 + 0A) 

0 


y 2 


0 

(-0.1220 + 0.4331A) 

(-0.0106 + 0A) 

(0.5745 - 1.8647A) 

0 


y 3 


0 

0 

0 

0 

(4- A). 


- y 4 - 


-0- 


(70) 


(71) 


(72) 


The new approximations for the lowest m = 2 eigenvalues from the Rayleigh quotients (eqs. (69)) are 


, = (V'^AV'!) = 0 
1 (Vi) T (BV' 1 ) 1.8647 


0.5745 


= 0.308079 


A^ = 


(V^) r (AV / 2 2 ) _ 1.30987 
(V^(BV^ 2 ) ~ 2.035 


= 0.64376 


(73a) 


(73b) 
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The first eigenvalue A^ is a very good approximation of the lowest eigenvalue for equations (39). An approximate 
check on Aj is possible since the exact value of A = Ai makes the determinant of the coefficients of equations (72) 
vanish. Substitution of A = Aj from equations (73a) into (72) makes the coefficients of the third row (or column) 
zero or near zero. An identically zero row (or column) is a sufficient condition for a determinant to be zero. 

The approximate check on the eigenvalue is the basis for the final step in the first iteration cycle for the 
second numerical example. The final step is contained in the next section of this report. 

The numerical example has n = 4 and m = 2, but the equivalence transformation for the example illustrates 
the matrix algebra for the general case. For the general case of m approximate eigenvectors of an n-dimensional 
problem, the matrix of the equivalent problem [Q r (A — AB)Q]Y = 0 differs from the original problem by only 
m rows and m columns. The m rows and columns are determ in ed by use of a permutation vector defined in the 
forward sweep procedure. In each of the m columns that are changed, (n — m) elements are simply the elements 
of (A — AB) times the corresponding approximate eigenvector. These elements are also stored for use as part 
of the Gram-Schmidt orthogonalization procedure. The coefficients in an m x m submatrix of the transformed 
problem are also found in the orthogonalization procedure. Therefore, starting with the eigenvector with the 
lowest subscript, the rules of matrix multiplication allow efficient computation and storage of the equivalence 
transformation and of the coefficients in the orthogonalization procedure simultaneously. 

Modified Sweep Procedures 

The modified sweep procedures are the last step in the iteration cycle for delayed division. The 
equivalence transformation and the modified sweep procedures form an inner iteration loop for improving 
the m approximate eigenvalues computed in the forward sweep and the m approximate eigenvectors computed 
in the back substitution. (See fig. 1.) 

After the equivalence transformation, the modified sweep procedures operate on the transformed problem 

(A' - AB')Y = 0 (74) 

where 

A' = Q r AQ 
B' = Q r BQ 

At this stage of the iteration, considerable information is available about the eigenvalue problem. The purpose 
of the modified sweep procedures is to use the information to reduce the number of computer operations in 
the iterative solution. Some information is related to exact properties of the problem and some information 
is related to only approximate properties. One exact property is that the transformed problem has the same 
eigenvalues as the original problem (A — AB)X = 0, whereas the eigenvectors are related through the change 
in variables X = QY. Because of the form of Q, only the elements of m rows and m columns of matrix A' 
differ from the corresponding elements of matrix A, with the remaining (n — m) 2 elements identical for A 1 and 
A. A similar statement holds for B / and B. 

Pertinent information on approximate properties is that the matrices A" defined by 

A" — A' - AfcB' (k = 1,2, ...,m) (75) 

contain at least one row and column that are nearly zero. The Gram-Schmidt orthogonalization procedure 
used in forming the matrix Q can be viewed as minimizing one column of each of the A" matrices in the 
least-squares sense. The nonzero elements in each of the m columns are a measure of the residual error in the 
orthogonalization . 

The matrices A" in equations (75) suggest that a shift in origin for each of the m eigenvalues would speed 
convergence of the iteration for solutions of equations (74). The modified sweep procedure makes the eigenvalue 
shifts simultaneously rather than sequentially. The simultaneous shift takes advantage of the fact that (n — m) 2 
elements of A and of B are unchanged by the transformation to A' and to B'. 

The operations for taking advantage of the unchanged elements are stored during the forward sweep 
procedure on the original problem (A — AB)X = 0. In matrix notation, the forward sweep procedure results in 

(L -1 + AM)(A — AB) = (R +• AS) + 0(A 2 ) (76) 
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The notation in equations (76) is an extension of the notation for Gaussian elimination subroutines. The 
matrix R is upper triangular and L is lower triangular. They are factors of A: 

A = LR (77) 

The matrix R is stored in the upper right-hand half of auxiliary matrix C during the forward sweep, whereas 
L -1 is stored in the lower triangular part of C as a product of elementary matrices. The matrices M and S 
in equations (76) are extensions of Gaussian elimination. Matrix S is upper triangular and is stored explicitly 
during the forward sweep, whereas M is lower triangular. Matrix M is stored as a product of elementary 
matrices derived by dropping quadratic terms in A at each elimination step in the forward sweep. 

After the equivalence transformation, it can be shown that ( n—m ) rows and columns of matrices (L -1 +AM) 
and (R + AS) are unchanged. Therefore, for small A, any eigenvalue shift can be applied to submatrices of 
(L _1 + AM) and (R + AS). If the shift is defined by 

A = A' + p (78) 

the approximation is 

(L- 1 + AM) = (L -1 + A'M) + pM (79a) 

(R + A'S) = (R + A'S) + pS (79b) 

The approximation in equations (79) is used in the modified sweep procedures for eliminating the first ( n — m ) 
variables in the transformed problem (eqs. (74)). Using the approximation in equations (79) avoids repeating 
numerical operations in the sweep for the transformed problem since the results of the operations have been 
stored during the forward sweep on the original problem. 

In the m columns of the matrices in the transformed problem (A' — AB')Y = 0 that are different from the 
original problem, the simultaneous shift is written explicitly as 

A = A fc + p fc (fc= l,2,...,m) (80) 

where A*, are current approximate eigenvalues stored after the equivalence transformation. The values of p& 
remain to be determined in the modified sweep procedure. In operating on any one of the m columns with 
elementary matrices of the form of equations (79a), the corresponding A^ value is substituted for A'. Because 
of the way elements are stored for the modified back substitution, a convenient identity (derived from eq. (80)) 
once a value is assigned to one of the p*. = p is 

Pi = {^k ~ ^i) 4" Pfc (81) 

The modified procedures can be thought of as partitioning the transformed problem. In an (n— to) x (n—m) 
submatrix of the transformed problem, convergence is quadratic in A^. In the rest of the transformed problem, 
the convergence is faster since it is quadratic in p^. 

The second numerical example (eqs. (39)) has n = 4 and m — 2. The value of n is too small to be a 
good illustration of equations (79), but the application of the simultaneous shift in equation (81) is clear. The 
transformed problem (A' — AB')Y = 0 is written out in equations (72) after the change in variables X = QY 
in equations (70). The eigenvalue shift in equation (80) is 

A — Ai + pi = 0.308079 4- pi (82a) 

A = A2 ~h p 2 — 0.64376 + P2 (82b) 

The m — 2 columns of the transformed problem that are shifted are columns 3 and 2. Equations (72) are 
rewritten 


(3 -A) 

(-0.4752 + 0.8369A) 
(-0.1220 + 0.4331A) 
0 


(0.06357 + 0.8369p 2 ) 
(0 - 2.035p 2 ) 

(-0.0106 + 0p 2 ) 

0 


(0.01143 + 0.4331pi) 

(-0.0106 + 0pi) 

(0- 1.8647pi) 

0 


0 ] [Yi] ro- 

0 Y 2 _ 0 

0 Y 3 ~ 0 

(4 — A) J LY4 J Lo_ 


(83) 
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The elimination of 
original problem (eqs. 
(see eqs. (40)) 


Yi and Y 4 in equations (83) is the same as the elimination of Xj and X 4 for the 
(39)). The forward sweep stores elements necessary to generate the elementary matrix 



r(l/3) + (A/9) 
0 
0 

. (0 + 0A) 


0 0 
1 0 
0 1 
0 0 


(0 + 0A) ' 

0 
0 

(1/4) + (A/16) J 


(84) 


The transformed equations are premultiplied by (C 14 r + AD 14 ). In forming elements of columns 1 and 4 of 
the product no shift in (C / 4 + AD^ 1 ) is used, whereas column 3 elements use the shift in equation (82a) and 
column 2 elements are computed from the shift in equation (82b). Which column to use in the shift with A 
is deduced from the permutation vector (eq. (53)). The result of the premultiplication of equations (83) by 
(C ^ 1 + AD/j 1 ), after dropping quadratic terms in A, /q, and /i 2 , is 


1 

(0.0257 + 0.3459/i 2 ) 

(0.0042 + 0.1604//!) 

O' 


'Yr 


O' 

(-0.4752 + 0.8369A) 

(0 - 2.035/i 2 ) 

(—0.0106 + 0 /q) 

0 


y 2 


0 

(-0.1220 + 0.4331A) 

(-0.0106 + 0 /i 2 ) 

(0 - 1. 8647/xi ) 

0 


y 3 


0 

0 

0 

0 

1 . 


.y 4 . 


. 0 . 


The elimination of Yi and Y 4 from equation (85) is completed with premultiplication by the elementary matrix 
(L^ + AMh): 


(l ^ 1 + am 14 ) 


1 0 0 
-(-0.4752 + 0.8369A) 1 0 

-(-0.1220 + 0.4331A) 0 1 
0 0 0 


0 

-(0 + 0A ) 
-(0 + 0 A ) 
1 


( 86 ) 


The premultiplication by + AM 14 J follows the same pattern as the premultiplication by (c i4 + AD, 4 j 

for eigenvalue shifts and truncation of quadratic terms. The result of the operations transforms equations (85) 
into 


'1 

(0.0257 + 0.3459/i 2 ) 

(0.0042 + O-1604/ii ) 

O' 


'Yf 


O' 

0 

(-0.0016 - 2.078/i 2 ) 

(-0.0097 + 0.0314/ii) 

0 


y 2 


0 

0 

(-0.0146 - 0.0653 /i 2 ) 

(-0.000048 - I. 868 / 11 ) 

0 


y 3 


0 

.0 

0 

0 

1 . 


_y 4 . 


. 0 . 


The modified sweep procedure ends at this stage after eliminating (n — m) = 2 variables. A total of m — 2 
equations are still coupled, but the coupling is weak. A relaxation type of algorithm is followed to compute 
the corrections m and /i 2 to the eigenvalues. Although the modified sweep has not produced a complete upper 
triangular matrix, the off-diagonal terms are small enough to neglect. Neglecting off-diagonal terms in the lower 
triangular matrix has the effect of reversing the order of the Gram-Schmidt orthogonalization procedure in the 
equivalence transformation. The orthogonalization procedure sweeps lower numbered approximate eigenvectors 
from higher numbered eigenvectors. The relaxation procedure sweeps higher numbered eigenvectors from lower 
modes. In the general case, the triangular equation form is in terms of equations reordered by the permutation 
vector. Equations (87) happen to be triangular without reordering the equations when the element in the 
third row and second column is neglected. The determinant of the triangular form vanishes when the diagonal 
elements are zero. The corrections to the eigenvalues have the following values that make the diagonal elements 
vanish: 


(-0.000047595) 

(-1.868297) 


= -2.5475 X 10 “ 5 


(-0.0016328) 

(-2.078) 


= -7.857083 X 10 -4 


( 88 a) 

(88b) 


After the corrections to the eigenvalues, /i/., are determined, the eigenvectors for the transformed problem are 
computed with the modified back substitution procedure. This modified procedure is similar to the original 


18 



back substitution procedure except that equation (81) must be used to relate the eigenvalue corrections that 
appear in m columns. 

The completion of the modified back substitution is the end of the first full iteration cycle for the second 
numerical example. The equivalence transformation and modified sweep procedures form an inner iteration 
cycle that can be repeated with stored data from the forward sweep procedure. 

The results for the second numerical example with one forward sweep and two modified sweeps are 
summarized in tables I and II. 

TABLE I. TWO LOWEST EIGENVALUES OF EQUATIONS (39) 


Operation 

Ai 

^2 

First forward sweep 

0.31101776 

0.68898224 

First Rayleigh quotient 

.30807949 

.64376187 

First modified sweep 

.30805401 

.64297616 

Second Rayleigh quotient 

.30797841 

.64310402 

Second modified sweep 

.30797867 

.64310399 

Exact solution 

0.30797853 

0.64310413 


TABLE II. TWO LOWEST EIGENVECTORS OF EQUATIONS (39) 



Value from — 


Component 

of 

First back 

First modified 

Second modified 

Exact 

eigenvector 

substitution 

back substitution 

back substitution 

solution 

un 

■ 


-0.44502 

-0.445042 



I 

1.0 

1.0 

UlZ 

-.822876 

■ 

-.801976 

-.801938 

u 14 


| ! 1 

0 

0 

U21 

-0.859602 

-0.826558 

-0.802697 

-0.801938 

u 22 

.548584 

.475492 

.445928 

.445042 

u 2?> 


1.0 


1.0 

^24 


0 


0 


The convergence of the method of delayed division is rapid. As might be expected from dropping quadratic 
terms in A at each step, the lowest eigenvalue and lowest mode shape converge faster than the second solution. 

Summary of Delayed Division 

Each step in the method of delayed division has been outlined separately in the previous sections of this 
report. In the actual computations the steps are independent, but certain data are stored for use in later 
steps. This section briefly reviews the data storage in the SPLIT subroutine as the storage relates to avoiding 
duplicate calculations. 

The storage requirements in the forward sweep subroutine parallel the storage in Gaussian elimination 
subroutines for the linear algebra problem AX = /. The main difference in storage requirements is that the 
forward sweep procedure operates on two matrices A and B in the eigenproblem (A — AB)X = 0. The matrices 
A and B are stored for later use in the equivalence transformation. The forward sweep procedure operates 
on two auxiliary matrices C and D. At the start of the forward sweep procedure, A = C and D = — B. As 
presented in equation (76), the forward sweep subroutine operates on matrices C and D to produce the form 

(L _1 + AM) (A - AB) = R + AS + 0(A 2 ) 

At the end of the forward sweep routine, elements of the upper triangular matrix R are stored in matrix C. 
Because of the pivoting strategy in the forward sweep, matrix R is not stored in the upper right-hand triangle of 
matrix C. A separate permutation vector is stored that can formally generate elementary permutation matrices 
that can form R as a matrix product. However, the permutation vector locates the elements of R so that R is 
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never arranged explicitly in triangular form. The elements of the elementary matrices whose product is L -1 
are also stored in C. Elements that are known to be one or zero are not stored. There is also room in matrix 
C to store factors in the determinant of A. Also, in the forward sweep, the elements of the triangular matrix 
S are stored in matrix D along with elements of elementary matrices necessary to form the lower triangular 
matrix M. 

The forward sweep is followed by the back substitution subroutine. The back substitution computes 
approximate eigenvectors by explicitly solving 

(R + A fc S)V* = 0 (90) 

The eigenvectors are stored in a matrix U with m columns and n rows. 

The back substitution subroutine is followed by the equivalence transformation. The transformed problem 
is Q t (A — AB)QY = 0. The matrix Q and its transpose Q 7 are not stored since the elements of Q are ones 
and zeroes and the m columns are stored in the matrix U. Most of the elements of the product Q r AQ are 
identical with matrix A. The m rows and columns that are different are stored in matrix C and similar terms 
from the product — Q r BQ are stored in matrix D. The m rows and columns are determined by the last m 
elements of the permutation vector. 

The choice of storage locations partitions the transformed problem Q 7 '(A — AB)QY = 0. A forward sweep 
eliminating the first (n — m) variables of the transformed problem would duplicate the same computations in 
the original forward sweep. However, the modified sweep procedure goes directly to the (n — m) rows and 
columns of matrices C and D and uses the (L _1 + AS) data that were stored on the forward sweep. Only 
m modal amplitudes remain to be eliminated on the modified sweep. The partitioning is such that quadratic 
terms in A^ are dropped in (n — m) columns of the modified forward sweep whereas quadratic terms in /if, are 
dropped in m columns, where is the correction on the current approximation for A*.. 

The modified back substitution computes approximate eigenvectors for the transformed problem. Only 
one vector W is stored at one time and is used to write the vector QWj.. The vector QW fc is stored in matrix 
U, updating the jth eigenvector Vfc for the original problem. Not storing all the transformed eigenvectors 
modifies Q so that the new approximate eigenvectors V*. are not computed exactly from the eigenvectors 
of the transformed problem. However, the approximation for the eigenvectors is consistent with the Gram- 
Schmidt orthogonalization procedure used to form the matrix U at the end of the forward sweep. The Gram- 
Schmidt procedure sweeps components of lower numbered approximate eigenvectors from higher numbered 
eigenvectors. The modified back substitution reverses the order and sweeps higher numbered eigenvectors from 
lower numbered approximate eigenvectors. 

The subroutines contained in the SPLIT subroutine make a complete iterative procedure. Convergence is 
problem dependent and is based on the assumption that m eigenvalues are small in absolute value compared 
with the remaining eigenvalues. In general, small eigenvalues converge faster than larger eigenvalues. The 
eigenvalues converge faster than the eigenvectors. The iteration is unaffected by multiple eigenvalues or closely 
spaced eigenvalues. 

These general conclusions are based on the theory and on results for some example problems presented in 
the next section of this report. After the example problems, the method of delayed division is compared with 
other methods for solving the symmetric eigenvector problem (A — AB)X = 0. 

Example Problems 

This section contains four eigenproblems solved with the method of delayed division. The problems are 
selected to illustrate different features of the method. The first problem is a column buckling problem modeled 
by finite differences. The column problem is a typical applied problem in which the matrix B of the problem 
(A — AB)X = 0 is not an identity matrix. 

The second example problem is a tridiagonal matrix problem posed by Wilkinson (ref. 4). Wilkinson uses 
the problem to show sources of numerical round-off errors in computing eigenvectors. However, the pivoting 
strategy in the method of delayed division is numerically stable and the computed eigenvectors are correct for 
the example problem. 

The third example problem requires more analysis. The problem is computing shear buckling loads for 
plates. The transverse plate deflection for a simply supported plate is represented in the numerical solution 
by a double Fourier sine series. The eigenvalue is proportional to the shear buckling load. The physics of the 
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problem dictates that for every positive eigenvalue there is a negative eigenvalue of the same absolute value. 
The occurrence of pairs of eigenvalues means that the strategy of dropping quadratic terms in delayed division 
fails for the shear problem. Two alternate methods for making linear terms dominant in the shear buckling 
problem are shown. 

The fourth example problem illustrates how an understanding of the method of delayed division allows the 
algorithm to be adapted to more sophisticated problems. The fourth problem applies the Lagrangian multiplier 
method to plate buckling. Through the use of Lagrangian multipliers, the order of the matrices A and B in 
the problem (A — AB)X = 0 is greatly reduced with no loss of accuracy in the final results. 


Column Buckling Problem 

The differential equation from linear theory for buckling of a column is 


El 


d 4 y p d 2 y 
dx 4 dx 2 


= 0 


(91) 


(See, for example, ref. 5.) The example problem is to compute the axial load P for a nontrivial solution of the 
differential equation subject to the boundary conditions for a clamped column. These boundary conditions at 
x — 0 and x — L are 

V = 0 



The numerical eigenvalue problem is formulated by approximating the differential equation by central 
differences at N interior nodes: 


Vi ~ 2 ~ ^Vi-1 + ®Vi ~ 4yi+l + Vi + 2 + KVi - 1 ~ 2 Vi + Vi+i) — 0 (i — 1, 2, ..., N ) 


The difference equation is written for N interior points. The eigenvalue A is 

PL 2 1 Ptt 2 

El {N + l) 2 ~ P e {N + l) 2 


(92) 


(93) 


where P E = ir^EI/L 2 is the Euler load for a simply supported column. The difference approximations for the 
boundary conditions are 

2/o — 0 and y N+1 = 0 
y- 1 = 2/1 and y^+2 = VN 

The unknowns on the left-hand sides of equations (94) are eliminated from equations (92), with N homogeneous 
equations left in N unknown at the interior nodes. 

The numerical results from the SPLIT subroutine are summarized in table III. The effect on the critical 
loads Pfc produced by varying N is summarized in table III. As N increases, the critical loads from the numerical 
eigenvalue solution approach the buckling loads from the exact solution of the differntial equation (91). The 
odd-numbered critical loads for fixed N can be checked from the exact formula 



Pk_ 

P E 


= 2 


^1 — cos 


fc7T + 7r\ / IV + 1 \ 2 

N+l) V tt ) 


(95) 


Equations (95) follow from an exact solution of the linear difference equations (92). The error tolerance on 
convergence for the iteration in the SPLIT subroutine was set to give six-digit accuracy for the critical loads. 
The lower buckling loads listed in table III hold this accuracy, and up to six eigenvalues were computed with 
good accuracy for each forward sweep. 
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TABLE III. CRITICAL LOADS AS A FUNCTION OF N 


(a) Odd-numbered loads 


N 

(Pi/PeY 

(Pi/Pe)' 

(Ps/PeY 

{P3/ P e)^ 

(Ps/PeY 

(p 5 /PeY 

4 

3.50056 

3.50056 





9 

3.87012 

3.87012 





10 

3.89242 

3.89242 

14.3338 

14.3339 



20 

3.97025 

3.97025 

15.5282 

15.5282 

33.6473 

33.6469 

50 

3.99494 

3.99494 

15.9192 

15.9192 

35.5928 

35.5921 

100 

3.99871 

3.99871 

15.9794 

15.9794 

35.8956 

35.8956 

oo 


4.00000 


16.0000 


36.0000 


*From SPLIT subroutine, 
f From equations (95). 


(b) Even-numbered loads (from SPLIT subroutine) 


N 

t P2/P E ) 

(Pa/Pe) 

(Pb/Pe) 

4 

5.874 



9 

7.544 



10 

7.652 



20 

8.0344 

23.0447 

43.89 

50 

8.1577 

23.9985 

47.44 

100 

8.1765 

24.1369 

47.99 

Exact 

8.1830 

24.2872 

48.19 


However, the example problem did reveal a limitation of the Gram-Schmidt orthogonalization procedure. 
The orthogonalization procedure is used to sweep lower modes out of approximate eigenvectors. The 
eigenvectors are columns in matrix Q in the equivalence transformation Q T (A — AB)Q. To compute 10 
or 20 eigenvectors in the matrix Q for the column buckling problem, the orthogonalization exhibited numerical 
instability. Although the eigenvalues are all distinct for the column buckling problem, the orthogonalization 
procedure failed to retain linear independence in the approximate eigenvectors. If linear independence is not 
maintained, the inner iteration loop fails to converge. 

The limitation imposed by the orthogonalization procedure on the number of eigenvalues and eigenvectors 
that can be computed with one forward sweep is not serious for postbuckling problems and for many other 
problems. In postbuckling problems, only a small number of nearly equal eigenvalues and eigenvectors are 
required for the analysis at any given load level. Therefore, the SPLIT subroutine is efficient for postbuckling 
problems. For problems in which more than a few eigenvectors are required, the orthogonalization procedure 
in SPLIT should be modified. The modification to the equivalence transformation is discussed in the section 
entitled Delayed Division and Other Methods. 

The efficiency of the iteration in SPLIT is indicated by computer times listed in table IV for executing 
different subroutines in the SPLIT subroutine. The times are CPU time in seconds on the CDC CYBER 
175 computer at the NASA Langley Research Center; however, the relative amount of time spent in each 
subroutine should be relevant for other computers. The forward sweep subroutine is executed once. It is 
the main subroutine in the method of delayed division. The forward sweep produces the factored form in 
equations (76) and includes the pivoting strategy and permutation vector described in the section entitled 
Forward Sweep. As might be expected from analogous results for Gaussian elimination subroutines on matrices 
of order n, the computer time for the forward sweep is roughly proportional to n 3 . 
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TABLE IV. APPROXIMATE COMPUTER TIME FOR SUBROUTINES ON CDC CYBER 175 


n* 

mt 

M t 

Subroutine CPU time for each call, sec, for — 

Forward 

sweep 

Equivalence 

transformation 

Modified 

sweep 

4 

2 

6 

0.020 

0.0075 

0.010 

9 

2 

7 

.025 

.0150 

.020 

10 

2 

5 

.040 

.0150 

.020 

10 

4 

9 

.040 

.0200 

.025 

20 

2 

5 

.170 

.0250 

.060 

20 ' 

4 

11 

.175 

.0450 

.110 

20 

6 

14 

.175 

.0700 

.150 

50 

2 

3 

2.255 

.1000 

.325 

50 

4 

5 

2.275 

.2100 

.645 

50 

10 

16 

2.285 

.5450 

1.525 

100 

2 

2 

17.451 

.3650 

1.280 

100 

10 

16 

17.856 

2.0400 

6.255 

100 

20 

16 

18.095 

4.3500 

12.300 


* Order of matrices A and B. 

^Number of computed eigenvalues and eigenvectors. 

■(-Number of iterations calling modified sweep subroutine and subroutine for orthogonalization and Rayleigh 
quotient. 


The equivalence transformation subroutine and the modified sweep subroutine are iterative and are each 
called by SPLIT a total of M times. The execution time for these iterative subroutines is approximately 
proportional to the product mn 2 , where m is the number of approximate eigenvectors in the equivalence 
transformation. 

The main conclusion that can be drawn from table IV is that the forward sweep procedure is efficient 
compared with the other subroutines. For matrices of the order of n — 100 or less, using to — 6 and then 
shifting to compute higher eigenvalues (to > 6) is an efficient strategy for SPLIT. 

Eigenvector Stability 

The column buckling problem revealed a numerical instability in the equivalence transformation if to, the 
number of computed eigenvectors per forward sweep, is too large. Wilkinson (ref. 4) discusses another type of 
numerical instability that can occur in computations for eigenvalue problems. This instability is manifested 
by errors in computed eigenvectors from known eigenvalues. Wilkinson recommends using a pivoting strategy 
for computation of eigenvectors. The complete pivoting in the forward sweep subroutine is numerically stable 
so that the eigenvectors computed by the SPLIT subroutine are numerically stable whenever the iteration for 
the associated eigenvalue converges. 

The example problem used here to illustrate the stable eigenvector computation is taken from reference 4. 


The problem has the simple tridiagonal form 

(aq — A)aq -F (3\X2 = 0 (96a) 

Pi*i-l + («* ~ A)x* + Pi+iXi+i = 0 (i = 2,3,4, ..., (n — 1)) (96b) 

Pn%n- 1 + fan ~ A)x n = 0 (96c) 

The particular tridiagonal form used in this example is 

Cii = —i+ 1 (*■— 1,2, .,.,«) (97a) 

Pi = 1 (97b) 


The terms a t in equations (97a) are shifted by 10 from Wilkinson’s problems so that the smallest Ai of 
0.74619418 for equations (96) with n = 21 corresponds to the largest Ai of 10.74619418 of reference 4. The 
exact eigenvectors are not affected by the shift. 
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Computed eigenvectors associated with Ai are listed in table V. The first column and the third column 
are almost identical. The first column is the first eigenvector from the SPLIT subroutine operating on the 
shifted problem [(A — Ail) — AI]X = 0, whereas the third column is reported by Wilkinson using pivoting 
for size on the matrix (A — Ail). The second column is from the SPLIT subroutine operating on the problem 
(A — AiI)X = 0. The second column is not quite as accurate as the first column, but the computed eigenvector 
is numerically stable. The fourth column is the eigenvector computed by Wilkinson directly from the matrix 
(A— Ail) without pivoting. The top components are correct, but round-off error makes the bottom components 
incorrect. 


TABLE V. STABILITY OF EIGENVECTOR COMPUTATION DUE TO 
PIVOTING IN DELAYED DIVISION 
[Ai = 0.74219418, first eigenvector u from delayed division] 



Uj for — 

i 

II 

O 

II 

l«< 

Pivoting* 

No pivoting* 

i 

1.0000000 x 10° 

1.0000000 x 10° 

1.00000000 

1.00000000 

2 

7.4619418 x 10“ 1 

7.4627843 x 10“ 1 

.74619418 

.74619420 

3 

3.0299994 x 10“ 1 

3.0284026 x 10“ 1 

.30299994 

.30299998 

4 

8.5902494 x 10“ 2 

8.5305094 x 10“ 2 

.08590249 

.08590260 

5 

1.8807481 x 10“ 2 

1.8413750 x 10“ 2 

.01880748 

.01880783 

6 

3.3614646 x 10“ 3 

3.2654781 x 10“ 3 

.00336146 

.00336303 

7 

5.0814709 x 10“ 4 

4.8526901 x 10“ 4 

.00050815 

.00051681 

8 

6.6594309 x 10“ 5 

6.2996079 x 10“ 5 

.00336146 

.00012346 

9 

7.7053623 x 10“ 6 

7.1409899 x 10“ 6 

.00000771 

.0004395 

10 

7.9828544 x 10“ 7 

7.3152313 x 10“ 7 

.00000080 

.003721 

11 

7.4882661 x 10“ 8 

6.6978590 x 10“ 8 

.00000007 

.03582 

12 

6.4181726 x 10“ 9 

5.6654776 x 10“ 9 

.00000001 

.3812 

13 

5.0644124 x 10“ 10 

4.3461741 x 10- 10 

Neglected 

4.442 x 10° 

14 

3.7025740 x 10“ 11 

3.1292977 x 10“ n 



5.624 x 10 1 

15 

2.5217695 x 10“ 12 

2.0630252 x 10“ 12 



7.686 x 10 2 

16 

1.6076278 x 10“ 13 

1.2922596 x 10“ 13 



1.128 x 10 4 

17 

9.6324699 x 10“ 15 

7.4583540 x 10“ 15 



1.1768 x 10 5 

18 

5.4443172 x 10“ 16 

4.1313141 x 10“ 16 



2.950 X 10 6 

19 

2.9121116 x 10“ 17 

2.1168481 x 10“ 17 



5.217 X 10 7 

20 

1.4783799 x 10“ 18 

1.0500629 x 10“ 18 



9.751 X 10 8 

21 

7.1260294 x 10“ 20 

4.8196742 x 10“ 2 ° 


■ 

1.920 x 10 10 


*From ref. 4. 


The fact that the example is a simple tridiagonal problem does not affect the general conclusion that 
the eigenvector computation in the method of delayed division is numerically stable. As with other iterative 
methods, the components of the eigenvectors computed by delayed division are not as accurate to as many 
digits as are the computed eigenvalues. 

In addition to the results already listed, results were computed for equations (96) with different input for 
m. As in the column buckling problem of the previous section, a value of to = 6 gives good numerical results. 
For to = 10, the orthogonalization procedure fails to preserve linear independence in the higher modes. 

Shear Buckling in Plates 

The third example is an applied problem on the shear buckling in orthotropic, simply supported plates. The 
matrix equations are generated from a Galerkin solution of the linear buckling equation with a double Fourier 
sine series used for the transverse deflection. The algebraic equations for the special case for an isotropic plate 
reduce to those solved by Stein and Neff (ref. 6) by power iteration. 

The linear buckling analysis is the first step in nonlinear postbuckling analysis. The equivalence transfor- 
mation applied in the eigenvalue computation is also useful for the nonlinear problem. However, this report 
only considers the linear buckling problem. 
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The forward sweep procedure in delayed division systematically retains linear terms in the eigenvalue while 
dropping quadratic terms. Delayed division fails when it is applied directly to the shear buckling problem 
because only quadratic terms appear in the forward sweep. The matrix D at the end of the forward sweep 
is zero. The analysis that follows shows that the terms for the shear load in the shear buckling problem are 
quadratic. The analytical results that shear buckling loads occur in pairs of opposite sign is not surprising on 
physical grounds, since a change in sign of a shear load merely indicates a change in direction. The insight 
provided by the analysis of the form of matrices A and B for the shear problem and by the fact that linear 
terms in the eigenvalue are retained in SPLIT allows modification of matrices A or B or both. 

The governing differential equation for the plate buckling problem is 


_ d^w , . d^w d^w 

Du a* + 2(Dl2 + 2Dm) feV + -^4 

d 2 w 

v l 4 > l\l 

y dy ‘ 2 


Ar d 2 w „ r d 2 w „ Ar 
— N x - - 0 + N v -^-^ + 2 N : 


dx 2 


xy 


dx dy 


(98) 


For the numerical example, all the coefficients in equation (98) are assumed constant. A Galerkin solution is 
sought in the form of a double Fourier sine series that satisfies the boundary conditions of simple support term 
by term. Thus, 


w 


. v'kx . airy 

— > > w v0 sm sin — — 

a b 

p=lq=l 


(99) 


The partial differential equation (98) is satisfied if the Fourier coefficients w vq are solutions of 


oo oo 
r= 1 S — 1 


IpqrsWrs — 0 


f P= 1,2,-,°°) 

\ (9 = 1 , 2 ,..., 00 ) 


( 100 ) 


where 


4 4 + 2 [(It) + ( w)l + qY (If)} - *** ~ K 


13 = 


K x = 


b 

N x b 2 


x ~ ix 2 D n 
N v b 2 


K - — 
Ky ~ it 2 D 


ll 


The coupling terms 'i vqr s are 


^ pqrs (p 2 _ r 2^q2 _ g 2 ) 

The eigenvalue to be determined is k : 


Ipqra = 0 (p = r or q = s) 
pqrs [1 - (— l) p+r ][l — (— l) 9+s ] 


k = 


4 

32 0K S 


TT“ 


(p ^ r or q ^ s) 


( 101 a) 

( 101 b) 

( 102 ) 


where K g = N X yb 2 /it 2 D\i. The symbol k is used for the eigenvalue because one matrix form of the problem 
defines X = k 2 for the eigenvalue. This matrix form follows from truncating the doubly infinite series in 
equations (99) and assigning single row subscripts to each pair of subscripts p and q. Thus, w pq = 27 when 
p — Pi and q = q- r Then, in the single subscript notation, equations (100) can be written in the standard form 
(A — kB)X = 0. The matrix A is a diagonal matrix with An — k pq for p — Pi and q = qn The matrix B 
contains the coupling terms j pqr3 defined in equations (101) and can be arranged in four submatrices. The four 
submatrices in matrix B are based on the form of the expression for 7 pqrs . Both (pi + rj) and ( qi + Sj ) must be 
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odd integers or B{j = — 7 vqrs = 0. The sum of two odd integers is even; therefore, the sum (p t + <&•) + (rj + Sj) 
is even or else B tJ = 0. For nonzero B^j, (p t + q^) and (ry + Sj) are both even or both odd. Therefore, the 
eigenproblem (A — /cB)X can be written in partitioned form as 


Aee 

0 

0 

0 1 


-* ee - 


' 0 

Boo 

0 

0 ' 


'X ee - 

0 

Aoo 

0 

0 


X 00 

= K 

B ee 

0 

0 

0 


V 00 

0 

0 

A-eo 

0 


*eo 

0 

0 

0 

Boe 


v eo 

0 

0 

0 

s 


.Xoe. 


0 

0 

B eo 

0 


Aoe. 


(103) 


where e and o indicate respectively even and odd values of the subscript integers. The partitioned problem 
decouples into two problems, 


Aee 

0 

and 

Aeo 

0 


0 ‘ 


A ee - 

= K 

' 0 

Boo 


X ee 

Aqo . 


LA 00 J 

. Bee 

0 J 


[v oo J 


0 1 


■x eo - 

= K 

' 0 

Boe 


'v eo - 

Aoe . 


Aoe. 

. B eo 

0 


.Xoe. 


(104a) 


(104b) 


The first problem (eqs. 104a)) contains eigenvectors with components w pq where {pi + qi) is even. The 
mode shape w from equations (99) for this case is symmetric with respect to 180° rotations about an 
axis perpendicular to the plate at its centroid. The second problem (eqs. (104b)) contains modes that are 
antisymmetric and (pj + q t ) is odd. 

The decoupled problems (eqs. (104)) are written as 

X ee = k {A~ 1 B 00 )X C0 (105a) 


A 00 X 00 — {BeeAgf} B 00 )X< 


- 1 : 


and 


X &0 — K,{A e J~ B 0 e )X 0 


(105b) 

(106a) 


A-oeAoe — K?(B eo A e J~ B 0 e)X 0e (106b) 

The method of delayed division for the SPLIT subroutine applied to equations (103) with A = k failed to 
converge because the final 2x2 matrix in the forward sweep procedure did not contain any A terms. However, 
the same subroutine with A — k, 2 gives good results when applied to equations (105b) and (106b). The 
results from SPLIT for the isotropic case verify the results obtained when Stein and Neff used power iteration 
(ref. 6 ). Stein and Neff retained 10 Fourier coefficients for each mode in computing the lowest eigenvalue for 
equations (104a) and (104b). The results from SPLIT retaining more Fourier coefficients do not change the 
buckling loads significantly. 

Another approach to solving equations (103) directly by delayed division is to use an eigenvalue shift. 
Let k = R + A, where ft is a nonzero approximate eigenvalue. Equations (103) are solved in the form 
[(A — kB) — AB]X = 0 with the SPLIT subroutine without any difficulty. 

Table VI contains a test case obtained from use of an eigenvalue shift on equations (103). The case is an 
isotropic plate with /? = a/b = 2.0. For this value of /3, the lowest buckling load associated with an asymmetric 
mode is nearly coincident with the lowest buckling load associated with a symmetric mode. Although the 
buckling loads are the same to an accuracy of two digits, there is no problem sorting out the two mode shapes 
when delayed division is used. 
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TABLE VI. FOURIER COEFFICIENTS FOR SHEAR BUCKLING PROBLEM 

[/? = a/b = 2.0] 


V 

g 

Lowest symmetric mode wpq from — 

Lowest asymmetric 
mode vjpq from 
SPLIT subroutine! 

SPLIT subroutine* 

Reference 7^ 

1 

l 

1.0000 

1.000 

0 

2 

2 

.3473 

.343 

0 

3 

l 

-.3373 

-.325 

0 

4 

4 

.0494 

.048 

0 

4 

2 

-.0482 

-.042 

0 

1 

3 

-.0318 

-.032 

0 

5 

1 

-.0080 

-.010 

0 

2 

4 

.0059 

.005 

0 

3 

5 

.0036 

.004 

0 

1 

5 

-.0019 

-.002 

0 

2 

1 

0 

0 

1.0000 

3 

2 

0 

0 

.3049 

1 

2 

0 

0 

-.2790 

4 

1 

0 

0 

-.1235 

2 

3 

0 

0 

-.0707 

4 

3 

0 

0 

.0318 

5 

2 

0 

0 

.0092 

3 

4 

0 

0 

.0044 


*K S = 6.563. 
tK s = 6.59. 
*K S = 6.610. 


The shear buckling problem illustrates a problem in which dropping nonlinear terms in the eigenvalue 
during the forward sweep can cause problems for the method of delayed division. Of the two approaches used 
to keep significant linear terms in the eigenvalue, the eigenvalue shift is the simplest to implement. The SPLIT 
subroutine only requires an input parameter to apply a shift. 

Lagrangian Multipliers 

The preceding example of shear buckling in plates illustrates the use of delayed division with an eigenvalue 
shift. The example in this section requires an eigenvalue shift and a new forward sweep with each iteration in 
the SPLIT subroutine. The eigenvalue shift with each iteration is essential because matrices A and B in the 
subroutine are not constant but are lambda matrices. A lambda matrix A is defined herein as a matrix whose 
elements are functions of a parameter A: 

Ai, = hM) ( 107 ) 

A typical problem containing lambda matrices requires calculation of characteristic roots A = X m for which 
the determinant of the matrix A vanishes. Characteristic roots of lambda matrices can be computed with 
delayed division. Fewer iterations are required for each root when delayed division is used compared with the 
number of iterations for determinant-search methods. The problem AX = 0 with A = fij(X) is rewritten for 
delayed division. The new problem is written with an iterative eigenvalue shift; convergence is achieved when 
the eigenvalue of the shifted problem approaches zero. 

The eigenvalue shift is 

^fc = ^fc-l + Mfc {k = 1,2, ...) (108a) 

where k is an iteration counter. The general form (A — AB)X = 0 for the method of delayed division is 


Aij — l) 


(108b) 




Bij — l) 




(108c) 
(108d) 

where /L denotes the derivative of /fy with respect to A. The smallest absolute value of the eigenvalue from 
the solution of (A — AB)X is taken as at each iteration step. 

A subroutine, LAMR, that solves equations (108) is described in reference 1 along with example problems 
with lambda matrices. The SPLIT subroutine is used herein to take advantage of symmetry (Ay = Aj, t ) in 
the example problem and to compute eigenvectors for multiple roots. 

The lambda matrix in the example problem is generated from the application of the method of Lagrangian 
multipliers to a plate buckling problem. The method of Lagrangian multipliers reduces a linear eigenproblem 
with many degrees of freedom to a nonlinear eigenproblem with only a few degrees of freedom. 

The particular problem treated herein is the buckling of an isotropic rectangular plate in compression and 
clamped on all four edges. An approximate solution for the clamped plate was obtained by Levy (ref. 8) in 
1942. Levy also reported results of earlier investigations on the problem. Thurston (ref. 9) worked the problem 
using Lagrangian multipliers as a special case of a sandwich plate with infinite shear stiffness. The buckling 
loads in references 8 and 9 were computed with a determinant search method. 

The analysis with Lagrangian multipliers starts when functions are found that satisfy the governing partial 
differential equation: 

“^ a2 “ " (109) 


V 4 tn + 


D dx 2 


0 


Consider functions of the form 


<t>pq = Qpq i x ) cos(qny/b) 


(P= 1) 2, ..., oo) 
{q = 0, 1, ..., oo) 


( 110 ) 


( 111 ) 


For each of the functions <j) pq which satisfy equation (109), the functions g = g vq {x) must satisfy the ordinary 
differential equation 

d 4 g a \2 d2 9 . t a \ 4 , ~N x a 2 d 2 g 

d?~ 2 {qfln) de +alqPn) +_ s _ 5? =0 

where ^ — x/a and /3 = a/b. The solution of equation (111) is 

g — ci cos(a + 6 )£ + C2 sin(a + <5)£ + C3 cos(a — 5)£ + 04 sin(a — <5)£ 
where a 2 — K pq (T[f3/2) 2 , 6 2 = a 2 — {q(3 tt) 2 , and 

—N x b 2 


KPC1 ~ 7T 2 D 


( 112 ) 


The solution g contains four constants of integration. The constants of integration are determined so that each 
function g pq {x) satisfies the following boundary conditions at x = 0 and x = a: 


9pq { x ) — 0 


(113a) 


dx 


= 0 


(113b) 


Since the boundary conditions (eqs. (113)) are homogeneous, the conditions determine a set of buckling 
coefficients K pq and associated mode shapes with one arbitrary constant of integration. 

After K pq and the modal function q pq {x) are determined for each value of p and q, the next step in the 
Lagrangian multiplier method is to assume the solution for the clamped-plate problem is 


w 


w pq g pq {x)cos{qTxy/b) 


(114) 


p= 1 q = 0 
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where w pq are constants to be determined. The solution w must satisfy the partial differential equation (109) 
plus the clamped-edge boundary conditions at x = 0 and x = a and at y — 0 and y = b of 


w = 0 


(115a) 

dw 

(x = 0 or x = a) 

(115b) 

at =0 

P=° 

dy 

(y = 0 or y = b) 

(115c) 


Since g pq (x ) satisfies equations (113), each term of the series for w satisfies all the boundary conditions at 
the edges x = 0 and x = a. The cosine terms in y are chosen to satisfy equation (115c). The two conditions of 
zero deflection at the unloaded edges y — 0 and y = b (eq. 115(a)) remain to be satisfied by the complete series 


X X w pq9pq( x ) = 0 (0 <x<a) 

p= 1 <j=0 
OO OO 

X X^p^wfaX -1 ) 9 = 0 (0<x<a) 

p= 1 q = 0 

or, through addition and subtraction of equations, the problem decouples into symmetric and antisymmetric 
solutions as follows: 


X X w pq9pq ( x ) =0 (0 < x < a) 

p= 1 <7=0,2 
OO OO 

X X w pq9pq( x ) =° (0 < X < a) 

P= 1 9=1,3 


(116a) 

(116b) 


The coefficients w vq must also be selected so that the series for w satisfies the partial differential equation (109) 
and gives 


V 74 , ( -N x \ d 2 w 7 r 2 ^ ^ d 2 g qwy 

v " + T £«>„(*. cos— = o 

p= 1 9=0 


(117) 


where K x = —N x b 2 /tt 2 D. The nondimensional buckling load coefficient K x is the parameter A in the lambda 
matrix for the clampled-plate problem. The lambda matrix is derived in the method of Lagrangian multipliers 
by use of a least-squares solution to reduce equations (116) and (117) to an algebraic problem. 

The functions g" q(x), the second derivatives of g p o(x), are a complete set of weighting functions for a 
least-squares solution of the boundary conditions (eqs. (116)). The boundary conditions become the algebraic 
constraint equations 


X X w PqI J ipq ~ 0 

P=1 9=0,2 

{i = 1 , 2 , ..., oo) 

(118a) 

OO OO 

^ 53 w P<l^ J ipq ~ 0 

p= 1 <2=1,3 

(z = 1 , 2 , ..., oo) 

(118b) 


The coefficients L^ pq in equations (118) are the definite integrals: 

ra ra 

Lipq = -k y o 9io( x ) 9pq{x ) dx = li g' lQ {x) g' pq {x) dx 
The constant l{ is a convenient normalizing factor: 


(119) 


5 = J>'^ 


dx 
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With the introduction of Lagrangian multipliers Hi and 7 t -, the algebraic problem is to minimize the function 
H: 


Jo Jo \ u ox j , =1 1 =Q 2 / . =1 \ 


H 


- N x d 2 


( 120 ) 


The function H must be minimized with respect to the coefficients w pq subject to the constraints in 
equations (118). 

The reason for the choice of (f) pq (eqs. (110)) as coordinate functions now becomes apparent. The quadratic 
form in the double integral simplifies because (j> pq axe orthonormal functions: 

rb ra ( at , 32 ,... \ 1 00 

Jo Jo w ^ w + J dxdy =2 ^^ w p^ k pi ~ K x ' ) ( 121 ) 

The orthogonality of the set of functions 


<Ppq = 9pqi x ) cos(qny/b) 


follows directly from the relations 


and 


f b 

/ cos(jny/b) cos(qiry/b) dy = 0 (j ± q) 

Jo 

ra ra ra 

J Q 9i q 9 pq dx = g'{ q g pq dx = - j g' iq g' pq dx = 0 (i ± p ) 


( 122 ) 

(123) 


The functions 4> pq are made orthornormal by proper choice of the arbitrary constant of integration in each 
function g pq {x). 

Substituting equations (121) into (120) and differentiating with respect to w pq gives the following set of 
equations: 

dH = 0 


dw. 


vq 


or 


{Kpq 

OO 

K x )w pq ^ ^ HjLjpq 0 
3 - 1 

/ (p = 1,2,. 
1 (? = 0, 2, ., 

8 8 

(124a) 

( K pq 

00 

— K x )w pq — ^ ^ ^ljLj pq = 0 

.7=1 

II II 

H— 1 

oo to 

8 8 

(124b) 


Solving equations (124) for w pq and eliminating w pq from the constraints (eqs. (118)) gives the following two 
decoupled sets of equations: 


OO OO T T 

^ipq^jpq 


E E E (K _ K ) 

j=l [p=lq=0,2 ^ Kpq 


E 

3=1 


OO OO T T 

^ipq^jpq 


E E 

P= 1 9=1,3 


(K pq - K x ) 


Hj = 0 (t = 1, 2, ..., oo) 

Ij = 0 (*' = l,2,„.,oo) 


(125a) 

(125b) 


The matrix of coefficients of the Lagrangian multipliers m and 7* in equations (125) is a lambda matrix with 
X — K x , the buckling load coefficient for the clamped plate. Each element of the lambda matrix is a double 
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infinite series rather than a polynomial in A as in the classic literature of lambda matrices. The coefficients of 
the multipliers m in equations (125a) are 


OO OO T T 

^ipq^jpq 


*y = E E 


p= l q= 0,2 


(K pq - K x ) 


(126a) 


A similar form appears in equations (125b). For nontrivial solutions of equations (125) the determinant of the 
lambda matrix must vanish. The iterative procedure in equations (108) is used in delayed division. Terms for 
Aij in the form of equations (126a) correspond to equations (108b). The corresponding B tJ terms in equations 
(108c) are 

oo oo r r 

d __T E ^ipq^jpq 

T-3 _ ir \2 

p= 1 <?=0,2 Kx) 

and A = AK X corresponds to the eigenvalue in equations (108d). 

For numerical solutions, the double infinite set of constraint equations (125) must be truncated. Also, 
the double series in equations (126) must be truncated. Theory shows (e.g., refs. 10 to 12) that by careful 
truncation either upper or lower bounds on K x can be computed numerically. Details of the theory are not 
explored herein except to note that the theory is useful in determining starting values for K x in the iteration 
using equations (108). 

The theory related to truncation provides guidance on the use of the SPLIT subroutine to compute the 
small eigenvalues A = A K x in the iteration sequence corresponding to equations (108). The forward sweep 
subroutine in SPLIT searches for large values of the coefficients Ay. However, K pq in the denominators of Ay 
and B t j in equations (126) leads to small values of Ay associated with even smaller values of By. To make 
the forward sweep subroutine eliminate higher modes first, the problem is “scaled.” The scaling is carried out 
formally with the change in variables 


(126b) 


x { = (-Bn) 1/2 Vi (i = l,2,...,n) 

so that (A* — AB*)Y = 0 is solved by the SPLIT subroutine for 

(127a) 

AL = (-%r 1/2 Ay ( -%)- i/2 

(127b) 

B*j = {-Bid-WBiji-Bjj)- 1 ' 2 

(127c) 


There is no danger of taking square roots of negative numbers since — B^ in equation (127b) is a sum of squares 
of real numbers. 

The matrix B* has negative ones for diagonal elements. The matrix A* has its largest absolute value 
elements on the diagonal, with small off-diagonal elements. The SPLIT subroutine picks out the smallest 
eigenvalue Ai = A K x from equations (127) and the iteration indicated by equations (108) is continued by 
shifting K x until Ai approaches zero. When Ai = 0 for equations (127), a critical value of K x is determined 
for a nontrivial solution of equations (125). The eigenvectors of equations (127) determine the Lagrangian 
multipliers and -yj, and equations (124) determine the coefficients w pq . The complete deflection mode shape 
is determined with equation (114). The number of degrees of freedom (i.e., the total number of w pq coefficients) 
can be an order of magnitude larger than the n constraint equations (125). The high accuracy of the final 
eigenvalues from a lambda matrix of relatively low order is the big advantage of the Lagrangian multiplier 
method. Table YII summarizes the results of Kx for the clamped plate with different number of degrees of 
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freedom and with different numbers of constraints. The iteration sequence of equations (108) converges 
rapidly for all cases. 


TABLE VII. BUCKLING COEFFICIENTS K x FOR CLAMPED PLATE 
(a) 20 constraints (eqs. 118); w vq = 1600 


/ 3 = a/b 

K x from — 

Method of 

Lagrangian multipliers 

Method of reference 8 

0.25 

64.0000 


.50 

19.3402 


.75 

11.6669 

11.569 

1.00 

10.0755 

” 10.074 

1.25 

9.2822 

9.25 

1.50 

8.3697 

8.33 

1.75 

8.1604 

8.11 

2.00 

7.9335 

7.88 

3.00 

7.4373 

7.37 

4.00 

7.3001 

7.23 


(b) f3 = a/b = 1.00; Lagrangian multiplier method 


Constraints 

Maximum 

V 

Maximum 

Q 

Total 

W PQ 

K x 

4 

400 

20 

20 

9.8670 

8 

600 

20 

30 

10.0211 

20 

1600 

40 

40 

10.0755 

20 

400 

20 

20 

10.1207 

20 

200 

20 

10 

10.1234 

20 

120 

20 

6 

10.1335 

20 

80 

20 

4 

10.1594 

8 

100 

10 

10 

10.1897 

8 

36 

6 

6 

10.2580 


The method of Lagrangian multipliers can be extended to more general problems. The solutions generated 
by program VIPASA for prismatic plate assemblies (refs. 13 and 14) satisfy periodic boundary conditions. The 
solutions with periodic boundary conditions obey orthogonality relations that give them a simple quadratic 
form similar to equations (121). Therefore, writing boundary conditions as constraint equations and using 
Lagrangian multipliers leads to lambda matrices similar to equations (125). The method of delayed division 
is a convenient tool for solving these equations for eigenvalues for the complete plate assembly under general 
boundary conditions. 

Delayed Division and Other Methods 

This report describes the method of delayed division for solving the symmetric eigenproblem (A— AB)X = 0. 
Other methods exist for solving the eigenproblem. The SPLIT subroutine that was written to implement 
delayed division also incorporates some operations from other methods. The purpose of this section is to 
note the operations used in the SPLIT subroutine that are common to other methods. After the common 
operations in the SPLIT subroutine are noted, a brief discussion follows on modifications that can be made 
to the SPLIT subroutine to make it more efficient for certain cases. One case in which the SPLIT subroutine 
should be modified is when the number of eigenvalues m computed per forward sweep is relatively large. 
Another separate case in which the SPLIT subroutine should be modified is when A and B are sparse, banded 
matrices. The novel features in the method of delayed division are in the algorithm for reducing the problem 
(A — AB)X = 0 to upper triangular form 

(L -1 + AM) (A - AB)X = (R + AS)X = 0 
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In the SPLIT subroutine, the reduction to the triangular form is contained in the forward sweep subroutine 
and the back substitution subroutine. The exact computation L -1 A = R allows storage of information on 
the factors of the determinant of matrix A. The factors can be used for Sturm sequence calculations when 
A and B are constant matrices and A or B is positive definite (ref. 15). The factors can also be used when 
the matrix A is a lambda matrix that obeys a Sturm sequence (refs. 2, 12, and 13). In addition to factoring 
matrix A, the method of delayed division retains linear terms in A. Therefore, because of the linear terms in 
A in the factorization, delayed division converges faster than determinant-search methods. As a faster solution 
than determinant search, the input for the SPLIT subroutine is m — 1 or 2 followed by an eigenvalue shift 
after each forward sweep. The solution of the plate buckling problem by use of Lagrangian multipliers is a 
good example of the advantages of using the SPLIT subroutine rather than a determinant search. 

The equivalence transformation X = QY in the SPLIT subroutine contains characteristics borrowed from 
inverse iteration. Bathe and Wilson review the literature and tabulate operation counts for inverse iteration 
in references 16 to 18. When m = 1, that is, when only one eigenvalue and corresponding eigenvector are 
computed, the relation between delayed division and inverse iteration is clearest. The single approximate 
eigenvector identified by the forward sweep and back substitution appears in matrix Q. The transformed 
matrix Q T (A — AB)Q contains the two terms in the Rayleigh quotient for the approximate eigenvalue. The 
additional feature with delayed division is that the pivoting in the forward sweep automatically computes 
the approximate eigenvector associated with the lowest eigenvalue. No additional logic is required to start 
the iteration. Since delayed division retains linear terms in A, the convergence of the inner loop containing 
the modified forward sweep and back substitution is faster than in simple inverse iteration. In the example 
problems of this report, acceptable engineering accuracy for the lowest eigenvalue and the associated eigenvector 
is obtained with one forward sweep and one modified sweep. 

When the lowest eigenvalues are closely spaced, simple inverse iteration converges slowly. Delayed division 
does not have any problem with several small eigenvalues. 

For several closely spaced, small eigenvalues, the version of inverse iteration that is appropriate is subspace 
iteration with Gram-Schmidt orthogonalization. The equivalence transformation X = QY in the SPLIT 
subroutine has much in common with this version of inverse iteration. The number of lowest eigenvalues in the 
SPLIT subroutine m corresponds to the dimension of the subspace. The forward sweep in SPLIT automatically 
selects the m approximate eigenvectors. The retention of linear terms in A speeds convergence. 

The Gram-Schmidt orthogonalization in the SPLIT subroutine is satisfactory when m is small, say on 
the order of six. For larger values of m, the assumption that A is small starts to break down in the sweep 
procedures and the Gram-Schmidt orthogonalization fails to preserve linear independence. The computer logic 
in the SPLIT subroutine needs to be modified after the equivalence transformation X = QY to treat the 
case of large values of m. The transformed matrix Q r (A — AB)Q contains an m x m submatrix that can 
be solved exactly for all its eigenvalues and its eigenvectors which will be orthogonal over the subspace. A 
branch in the SPLIT subroutine with an exact solution over the subspace would be analogous to the method 
of inverse iteration with calculation of operator projections. Here again, as in the case for computing only one 
eigenvector, the logic in the forward sweep of the delayed division automatically computes the m vectors that 
span the subspace and speeds convergence. 

The partitioning of the matrix Q r (A— AB)Q involved in the exact m dimensional subspace solution also has 
features similar to the Guyan reduction (ref. 19). Again the pivoting strategy in the forward sweep automates 
the reduction. Implicit in the Guyan reduction is the assumption that matrix A is positive definite. When A 
is positive definite, all the factors in its determinant are positive so that it is obvious that the pivoting in the 
forward sweep drives small eigenvalues of matrix A into the subspace. The retention of linear terms in A in 
the forward sweep is a practical method of implementing the algorithm suggested by Kidder for modifying the 
Guyan reduction (ref. 20). 

The partitioning suggested by the m dimensional subspace solution has already been programmed in the 
SPLIT subroutine. A branch solving the m dimensional subspace problem exactly should be simple to add. 

The case for which A and B are large, sparse matrices requires a completely new subroutine in order to 
apply the method of delayed division. Pivoting is required to maintain numerical stability. Where A is positive 
definite, pivoting on the main diagonal is stable, does not increase the bandwidth of the reduced equations, 
and retains symmetry. Details of the pivoting strategy and storage requirements are beyond the scope of this 
report; however, the approach suggested by Gupta (ref. 15) for the storage problem is feasible for the algorithm 
in delayed division. 
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This discussion has shown that combining the advantages of the method of delayed division with other 
methods may be useful. As noted by many authors, the best combination of methods is problem dependent 
and there is no “best” solution to the eigenproblem (A — AB)X = 0. 

Concluding Remarks 

The method of delayed division for efficiently solving the symmetric eigenproblem for small eigenvalues has 
been developed. Delayed division converges rapidly when the number of computed eigenvalues is small. The 
method computes the associated eigenvectors of the problem directly. The computation of the eigenvectors is 
numerically stable and has been programmed in a subroutine called SPLIT. 

The SPLIT subroutine has been programmed to apply delayed division when the number of computed 
eigenvalues is small and A and B are dense matrices. The subroutine could be revised for cases with a larger 
number of computed eigenvalues. A new pivoting strategy is required before the method is efficient for large, 
sparse matrices. 


NASA Langley Research Center 
Hampton, VA 23665-5225 
October 7, 1985 
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Symbols 


P 

constant 

A,B,C,D,G, 

L,M,R,S 

A*,B* 

a 

Gj, 1 . ^5 *1 , 

coefficient matrices 

scaled coefficient matrices 

constant; plate length 

elements in zth row and yth column 

Q 

q 

u 

u 

matrix in eigenvalue transformation 
constant 

matrix of eigenvectors 
component of U 

l J ’ C J 1 
c ij i , 

of coefficient matrices 

V 

matrix of approximate eigenvectors 

ffij 


V 

component of V 

b 

constant; plate width 

w 

transverse deflection of plate 

c 

constant 

X 

vector 

D pq 

function in series solution for w in 
plate problem 

X 

axial coordinate of column; 
Cartesian coordinate 

El 

bending stiffness of column 

Y 

vector in equivalence transformation 

9pq 

function of x in assumed form of 
solution for plate problem 

y 

transverse deflection of column; 
Cartesian coordinate 

9pq 

derivative of g pq with respect to £ 
in plate problem 

*-b' > Pi 

nonzero elements of tridiagonal 
matrix 

H 

function to be minimized in La- 
grangian multiplier method for 

a, /3,8 

constants in plate problem 


plate problem 

A 

determinant of 2 X 2 matrix 

I 

identity matrix 

V 

biharmonic operator 

K pq 

buckling coefficient for plates in 

K 

eigenvalue in plate problem 


compression 

A 

diagonal matrix containing 

Ex, Ky 

buckling coefficients for clamped 


eigenvalues 


plate 

A 

eigenvalue 

L 

length of column 

P 

shifted eigenvalue 

Lipq 

coefficient in plate problem 

Pj, 1 j 

Lagrange multipliers in plate 

m 

number of computed eigenvalues 


problem 


and corresponding eigenvectors 

z 

= x/a, nondimensional variable in 

N 

number of interior stations for 


plate problem 


column finite-difference equations 

4>pq 

function in series solution for w in 

N x 

stress resultant in plate problem 

plate problem 

n 

total number of rows and columns 

Subscripts: 



in matrices A and B 

I,J 

constants in elimination of X 

0( A 2 ) 

order of A 2 

i,3,k,p, 

constants 

P 

axial load on column 

q,r,s 


P 

permutation matrix 

Primes indicate matrices or elements after ele- 

P E 

Euler buckling load for simply 
supported column 

mentary operations in analysis or in computations. 
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